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ABSTRACT 


The  performance  of  Loran-C  receivers  is  severely  limited  by 
contamination  introduced  by  multipath  propagation  effects.  Echo  removal 
by  discrete  generalized  linear  filtering  is  considered  as  a  method  of 
reducing  tlie  effects  of  the  contamination.  The  signals  can  be  described 
as  consisting  of  repeated,  knovm,  composite  signals  with  short  echo 
arrival  times,  in  noise.  Since  the  signals  are  repeated,  and  stable, 
they  can  be  averaged  to  increase  the  effective  signal  to  noise  ratio. 
The  homomorphic  deconvolution  technique  is  applied  to  simulated  signals 
and  a  determination  of  the  averaging  time  needed  for  satisfactory 
performance  is  made. 

Aside  from  the  Loran-G  application,  the  report  presents  a  detailed 
study  of  the  computational  aspects  involved  in  deriving  the  complex 
cepstrum.  Since  the  signal  is  known,  the  effects  of  noise  and  aliasing 
can  be  precisely  identified.  In  particular,  conventional  phase 
um-frapping  methods  are  questioned  for  certain  classes  of  signals. 
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Chapter  1 
INTRODUCTION 

The  digital  signal  processing  technique  referred  to  a  Homomorphic 
Deconvolution  has  met  with  some  success,  in  a  number  a  diverse  appli- 
cations, in  suppressing  signal  contamination  introduced  by  a  multi- 
path  propagation  medium.  Another  possible  application  is  to  Loran-G, 
a  low  frequency  electronic  navigation  system,  in  which  the  performance 
is  severely  limited  by  the  effects  of  such  contamination.   It  is  the 
purpose  of  this  thesis  to  study  the  applicability  of  the  technique  to 
various  aspects  of  Loran-C  signal  processing. 

In  this  Chapter,  a  brief  description  of  the  Loran-C  system  and 
Homomorphic  Deconvolution  is  presented.   In  Chapter  2,  the  task  of  the 
homomorphic  filter  in  the  Loran-C  case  is  foirmulated.  In  Chapter  3> 
the  filtering  technique  is  developed  ajid  the  results  of  simulated 
filtering  in  the  noiseless  case  aire  presented.  Chapter  ^■   presents  the 
theory  and  a  simulated  example  of  the  effects  of  signal  distortion. 
Chapter  5  includes  a  theory  of  the  effect  of  white  gaussian  noise  on 
the  processing  and  Includes  examples  to  verify  the  theory.   In 
Chapter  6,  the  effects  of  the  various  assumptions  and  simplifications 
made  are  discussed  and  the  results  of  the  preceeding  chapters  are 
examined  to  determine  the  applicability  of  the  method.  The  examples 
presented  are  simulations  performed  on  an  IBM  370  general  purpose 
computer. 
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1.1  Loran-G  System  Description 

Loran-G  is  an  electronic  navigation  system  which  belongs  to  a 
broad  class  of  systems  referred  to  as  hyperbolic,  a  characterization 
based  on  straightforward  geometric  considerations.   If  radio  signals 
aire  simultaneously  transmitted  from  two  different  locations,  they  will, 
in  general,  arrive  at  some  arbitrary  third  location  at  two  different 
times.  To  the  extent  that  the  speeds  of  propagation  along  the  two 
paths  are  equal  and  known,  a  measurement  of  the  difference  in  times 
of  arrival  provides  a  corresponding  measure  of  the  difference  in 
distances  involved.  The  locus  of  points  with  a  constant  difference 
in  distance  from  two  reference  points  is  a  hyperbola.  A  receiver 
which  determines  this  time  difference  and  identifies  the  proper 
reference  points  provides  information  sirfficient  to  define  a  navi- 
gational line  of  position.  Two  or  more  such  lines  of  position  can 
be  used  to  obtain  a  fix  -  the  location  of  the  receiver, 

Vaxious  hyperbolic  systems  with  origins  in  the  late  1930 's  have 
evolved  into  today's  alternatives.  Of  these,  it  may  be  said,  the  one 
which  meets  the  most  general  operational  requirements  is  the  Loran-G 
system.  The  advantages  of  Loran-G  can,  in  short,  be  attributed  to 
the  following, 

1.  Use  of  a  frequency  band  (90  to  110  kHz)  which  provides 
for  extended  range  and  exhibits  excellent  propagation 
properties, 

2,  Use  of  a  pulsed  signal  format.  This  increases  average 
transmitted  power  and  permits  signal  coding  to  facilitate 
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transmitting  station  identification  and  to  reduce  the  effects 
of  interference, 
3.  Improved  measurement  precision  obtained  by  coarse  measure- 
ment of  the  envelope  followed  by  fine  tune  measurements  based 
on  the  signal  rf  carrier. 
In  Loran-G,  the  hyperbolic  navigation  concept  is  realized  in  a 
somewhat  complex  manner.  There  are  a  number  of  modifications  to  the 
simple  scheme  as  previously  described.  While  the  modifications  maJce  the 
system  appear  involved,  they  are  seen  to  be  well  founded  and  straight- 
forward when  systematically  examined, 

A  first  consideration  is  the  operation  of  the  transmitting  stations. 
These  are  grouped  into  chains  whose  configurations,  subject  to  practical 
constraints,  provide  geometrically  optimal  accuracy  and  coverage.  An 
example  of  a  simple  chain  to  analyze  is  the  so-called  triad  which 
consists  of  stations  denoted  master,  slave  X  and  slave  Y,  The  trans- 
missions consist  of  groups  of  pulses  on  a  100  kHz  carrier  frequency: 
the  master  transmits  nine  pulses  in  a  group  while  the  slaves  each 
transmit  eight  pulses  in  a  group.  In  practice,  the  stations  do  not 
transmit  simultaneously.  The  slave  stations  incorporate  known  trans- 
mission delays  relative  to  the  master  so  that,  within  the  service 
area  of  a  triad,  reception  is  typically  as  depicted  in  figure  1-1, 
Transmission  timing  is  controlled  so  that  the  groups  are  received  in  the 
order  and  intervals  shown.  The  transmission  delays  guarantee  that  there 
will  be  no  overlap  of  the  groups.  In  this  example,  the  transmission 
sequence  is  repeated  every  50 i 000  microseconds  -  the  chain  pulse 
group  repitition  interval.  Different  chains  are  distinguished  by 
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Figure  1-1  Loran-G  Chain  Pulse  Transmission  Scheme 


different  group  repitition  intervals.  Consequently,  by  means  of  a 
correlation  type  operation  applied  over  a  number  of  repitition  intervals, 
a  receiver  can  discriminate  between  signals  transmitted  from  adjacent 
chains. 

The  stations  also  employ  a  pulse  phase  code  in  the  transmission. 
This  allows  enhanced  discrimination  between  signals  from  competing 
chains  and  overcomes  some  of  the  problems  introduced  by  skywave 
contcimination. 

The  term  skywave  refers  to  the  multipath  propagation  phenomenon 
depicted  in  figure  1-2.  The  high  accuracy  and  coverage  advantages  of 
Loran-C  depend  on  the  use  of  measurements  obtained  from  signals 
arriving  via  the  stable  groundwave  propagation  mode  at  extended  ranges. 
The  reception  of  skywaves,  whose  arrival  times  axe  highly  dependent 
upon  unpredictable  ionospheric  properties,  causes  severe  signal 
processing  problems.  The  nature  of  the  important  properties  of  the 
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Figure  1-2  Loran  Propagation  Paths 


multipath  medium,  at  Loran-G  frequencies,  is  illustrated  in  figures 
1-3  and  1-4.  It  can  be  seen  that  the  times  of  arrival  can  be  as  short 
as  about  35  microseconds  after  the  arrival  of  the  groundwave  for  strong 
first  hop  skywaves  at  extended  ranges.  Additionally,  arrival  times  can 
be  delayed  as  long  as  about  I5OO  microseconds  for  higher  order  hops 
vmder  some  conditions.  An  explanation  of  how  these  effects  are 
overcome  requires  a  further  description,  of  the  transmitted  signal 
format  and  a  brief  description  of  receiver  detection  methods. 

All  transmitting  stations  attempt  to  generate  pulses  of  the 

2  -at 
same  form  which  can  be  described  as  t  e    with  a  such  that  a  peak 

occurs  65  microseconds  after  the  beginning  of  the  pulseo  Such  a  pulse 

has  an  effective  duration  of  250  to  300  microseconds.  An  example 

which  combines  the  two  most  deleterious  effects  of  the  multipath 

contamination  is  illustrated  in  figure  1-5,  For  clarity,  the  magnitudes 
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First  Pulse  Skywaves 


Second  Pulse  Skywaves 


^First  Pulse  Groundwave 
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Figure  1-5  Composite  Received  Signal 


of  the  separate  components  of  the  signal  have  been  superimposed. 

The  effects  of  the  contamination  can  be  grouped  into  two 
categories.  The  groundwaves  are  overlapped  by  skywaves  associated  with 
the  same  transmitted  pulse.   These  skywaves  correspond  to  low  order 
hops  and  may  be  many  times  larger  than  the  groundwave.  The  groiindwaves 
are  also  overlapped  by  much  delayed,  high  order  hops  associated  with 
the  preceeding  transmitted  pulse.  These  skywaves  will  be  fairly 
well  attenuated  and  smaller  than  the  groundwaves  they  overlap. 

The  problems  of  long  delayed  skyi^ave  interference  is  easily 
overcome  by  use  of  a  pulse  phase  coding.  The  phase  of  the  trans- 
mitted signal  carrier  is  systematically  varied  in  180  steps  from 
pulse  to  pulse  in  accordance  with  the  scheme  shown  in  figure  1-6, 
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First 
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Third 

same  as  first 
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Fourth 

same  as  second  . 
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etCc 

etc.           1          etc.       1 

Figure  1-6  Pulse  Phase  Code 

As  an  example  of  the  use  of  the  code,  consider  the  effect  of  the 
sixth  hop  skywave  in  figiire  1-5  which  will  cause  problems  hy  over- 
lapping the  next  groundwave.   It  can  be  verified  from  figure  1-6  that 
the  code  can  be  represented  as  a  signal,  s(t),  vfith  the  property  that 
s(t)  and  s(t-T)  are  orthogonal  vrhen  averaged  over  two  pulse  group 
repitition  intervals,  where  T  is  a  multiple  of  1000  microseconds. 
Since  the  effect  of  the  sixth  hop  skywave  is  to  create  a  signal 
s(t)  +  A  s(t-T),  if  the  received  signal  is  correlated  with  s(t),  this 
type  of  interference  is  removed. 

The  effects  of  the  interference  caused  by  low  order  skywaves 
are  not  so  easily  overcome.  Use  is  made  of  the  fact  that  the 
groundwave  travels  the  most  direct  path  and  is  therefore  guaranteed 
to  arrive  before  any  of  the  skywaves  associated  with  the  same  trans- 
mitted pulse 0  Consequently  the  first  portion  of  the  groundwave  is 
free  of  contamination  of  this  type  and  can  be  used  to  make  the  time 
measurements.  This  portion,  however,  can  not  be  assumed,  a  priori, 
to  be  much  more  than  the  first  30  microseconds  of  the  pulse. 

The  measurement  technique  proceeds  in  two  steps.  First  it  is 
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Figure  1-7  Envelope  Deriver  and  RF  Samplers 
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attempted  to  identify  some  point  on  the  leading  edge  of  the  pulse. 
Since  this  point  must  be  within  the  first  30  microseconds  after  the 
beginning  of  the  pulse,  it  is  clear  that  this  is  a  severe  constraint 
on  a  pulse  which  has  an  effective  bandwidth  of  only  10  kHz. 

A  typical  technique  employed  is  to  process  the  pulse  envelope 
so  as  to  produce  a  zero  crossing  at  the  25  microsecond  point.   This 
is  accomplished  by  delaying  the  received  signal  5  microseconds  and 
scaling  it  to  produce  what  is  represented  by  the  dashed  line  in 
figure  l-7(a).  This  is  then  subtracted  from  the  original  signal  - 
the  solid  line  in  figure  l-7(a)  to  produce  the  derived  envelope 
of  figure  l-7(b).  The  zero  crossing  is  then  detected  and  sampling 
strobes  are  centered  at  the  estimated  25  microsecond  point  of  the 
rf  signal  as  shown  in  figure  l-7(c). 

With  the  strobes  located  as  shown,  the  middle  one  will  obtain 
a  zero  sample  value  while  the  other  two  will  obtain  values  with  a 
known  ratio.   If  the  strobes  were  not  centered  properly,  different 
values  would  be  obtained,  producing  an  error  signal  which  can  be 
used  to  reposition  the  strobes, 

In  practice,  of  course,  the  signals  are  noisy  so  that  both 
steps  in  the  procedure  must  be  accomplished  via  measurements  on  a 
number  of  pulses  -  as  many  as  6^  or  128  or  more  -  depending  on  the 
signal  to  noise  ratio.  Typically,  the  most  difficult  step  in  the 
procedure  is  the  envelope  zero  crossing  deriving o  The  zero  crossing 
measurement  must  be  accurate  to  within  5  microseconds  or  the  fine 
tune  measurement  that  follows  will  be  made  on  the  wrong  cycle  of 
the  caxrier.  In  this  regard,  a  final  chairacteristic  of  the  signals 
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should  be  presented. 

The  propagation  paths  involved  do  constitute  dispersive  media 
so  that  the  signals  received  are  not  quite  identical  to  those  trans- 
mitted. The  distortion  in  the  shape  of  the  envelope  is  slight  -  so 
much  so  that  even  though  there  are  different  propagation  paths 
involved,  the  skywaves  may  be  modelled  as  perfect  echoes  of  the 
groundwave.  The  major  effect  of  the  propagation  paths  is  the  intro- 
duction of  a  phase  shift  of  the  pulse  envelope  relative  to  the  carrier. 

In  effect,  the  envelope  and  carrier  propagate  at  different  velocities. 
Since  the  final  time  difference  measurement  is  based  on  the  carrier 
times,  the  carrier  propagation  times  are  calibrated  for  measuring 
distances.  While  the  phase  shift,  for  the  groundwave,  is  typically 
not  enough,  by  itself,  to  cause  a  cycle  ambiguity  in  the  detection, 
it  does  require  that  the  envelope  deriver  have  accuracy  somewhat 
better  than  5  microseconds. 

With  the  preceeding  as  background,  some  aspects  of  system  per- 
formance can  be  presented.  Inaccuracies  can  be  introduced  by  any 
number  of  sources  in  the  system.  These  can  be  loosely  grouped  into 
a  few  categories  and  statistically  described.  Receiver  error  will 
have  a  standard  deviation  of  from  20  to  100  nanoseconds  depending  on 
the  degree  of  sophistication.  Chain  synchronization  includes  the 
results  of  many  effects  and  may  be  described  as  introducing  an 
error  with  a  standard  deviation  of  about  50  nanoseconds.  Noise 
induced  errors  can  have  a  standard  deviation  of  from  about  20  nano- 
seconds for  0  db  SNR  to  about  250  nanoseconds  for  -20  db  SNR,  The 
implications  of  these  timas  can  be  seen  from  an  example.   If  a  total 
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error  of  about  100  nanoseconds  is  present,  this  corresponds,  via  a 
speed  of  propagation  conversion, to  a  distance  error  of  about  30 
meters. 

The  numbers  presented  above  for  the  error  statistics  were 
purposely  given  without  mention  of  important  descriptive  j^arameters. 
For  example,  a  figure  for  the  error  caused  by  atmospheric  noise  is 
completely  meaningless  without  a  detailed  description  of  the 
particular  receiver  detection  procedure  -  particularly  the  sample 
averaging  or  integration  time  involved  (50  seconds  in  the  case  of  the 
above  figures).  The  geometry  of  the  chain  and  the  receiver  position 
must  also  be  specified:  the  fix  obtained  from  lines  of  position  inter- 
secting at  right  angles  can  be  expected  to  have  an  error  almost  an 
order  of  magnitude  smaller  than  if  the  lines  were  nearly  asymptotic 
and  at  extreme  ranges.  Without  going  into  great  detail,  the  per- 
formance can  be  simply  summarized, 

1,  Position  errors  vary  from  a  few  tens  of  meters  to  not 
more  than  about  200  meters  within  the  coverage  region, 

2,  The  only  hyperbolic  navigation  systems  that  are  competi- 
tive in  performance  have  a  coverage  range  limit  of  less 
than  200  miles,   Loran-G  has  a  nominal  coverage  range 

of  about  1000  miles, 

3,  The  receiver  generally  contributes  the  smallest 
percentage  of  the  total  error. 

In  spite  of  the  third  property,  improved  receiver  performance  is 
constantly  strived  for  since  it  indirectly  effects  many  of  the 
other  sources  of  error.  Chain  synchronization  is  a  single  example  of 


-22- 

such  a  source.  Transmission  timing  is  controlled  by  passive 
Cesium  Beam  frequency  standards  located  at  each  station.  As  precise 
as  these  may  be,  there  will  be,  over  a  period  of  time,  relative 
phase  drifts  between  the  stations  of  a  chain.  These  phase  drifts 
are  observed  by  monitoring  stations  which  use  the  most  sophisticated 
receivers  and  periodic  corrections  are  ordered  to  maintain  syn- 
chronization to  within  100  nanoseconds.  Clearly,  if  more  sensitive 
receivers  were  available,  this  monitoring  process  could  be  made 
more  responsive  to  oscillator  phase  discrepancies. 

Speculation  on  means  of  achieving  improved  processing  schemes 
is  tempting  and  many  faceted.  One  aspect  to  be  considered  is  that 
an  optimum  receiver  need  not  totally  restrict  the  measurements  to  be 
made  in  the  portion  of  the  ground wave  that  is  perfectly  free  from 
skywave  contamination.  If,  for  example,  the  sampling  strobes  are 
centered  at  the  30  microsecond  point  rather  than  at  the  25  micro- 
second point,  the  signal  to  atmospheric  noise  ratio  would  be  im- 
proved. In  such  a  case,  it  might  result  that  the  third  strobe 
would  be  sampling  a  signal  contaminated  by  a  skywave.   In  many  cases 
however,  the  skywave  would  be  of  such  a  low  power  level  at  the 
sampling  point  that  the  increase  in  total  noise  power  -  atmospheric 
noise  plus  skyi^ave  interference  -  is  more  than  offset  by  the  in- 
crease in  signal  power.  The  implementation  of  the  resulting 
adaptive  time  filter  would  require  tlie  knowledge  of  skywave  relative 
magnitude  and  arrival  time  so  that  a  decision  could  be  made  whether 
or  not  to  re-position  the  sampling  points, 

A  less  optimal,  but  still  improved,  receiver  would  make 
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measurements  on  the  skywave  free  portion  of  the  signal  -  but  at 
later  times  in  the  many  cases  in  which  skyviave  arrival  times  are 
much  more  than  30  microseconds  after  the  beginning  of  the  groundvrave. 
In  this  case,  the  receiver  requires  only  the  knowledge  of  the 
skywave  arrival  time.  More  ambitious  attempts  have  been  made  to 
implement  a  total  power  receiver.   In  this  receiver,  an  attempt  is  made 
to  digitally  implement  a  scheme  to  measure  skywave  arrival  times, 
magnitudes  and  phases,  and  then  to  subtract  them  from  the  received 
signal  thereby  greatly  Increasing  the  signal  to  noise  ratio. 

Another  improvement  that  can  be  thought  of  would  be  allowed  by 
a  receiver  capable  of  making  precise  measurements  on  the  envelope  of 
the  groundwave.  This  would  overcome  the  problems  introduced  by  the 
different  propagation  speeds  which  might  cause  cycle  ambiguity  in  the 
measurement.  Even  if  a  measurement  based  on  the  envelope  might 
not  have  quite  the  precision  of  the  rf  measurement  and  therefore 
not  be  of  use  in  a  navigation  receiver,  it  would  be  useful  in 
calibrating  regions  in  which  the  envelope  -  carrier  phase  shift 
is  significant. 

The  investigation  of  an  alternate  method  of  accomplishing  some 
of  these  goals  is  the  subject  of  this  report.  The  method  -  homo- 
raorphic  deconvolution  -  will  now  be  presented, 

1.2  Honomorphic  Deconvolution 

In  1965»  A.V.  Oppenheim  presented  a  theory  of  generalized  super- 
position which  could  be  applied  to  certain  classes  of  non-linear 

2 

systems.   A  simplified  description  of  the  theory  can  be  given  via 
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an  example.  Consider  processing  a  signal  consisting  of  tvjo  or  more 
components  which  have  been  combined  in  some  non-linear  fashion.  By 
means  of  certain  elementary  operations,  it  may  be  possible  to  transform 
the  non-linear  combining  operation  into  a  linear  one.   If  the  components, 
when  so  transformed,  show  a  certain  distinctiveness,  it  may  be  possible 
to  process  them  by  conventional  linear  techniques.   One  of  the  components 
can  be  enhanced  relative  to  the  others,  or  the  components  can  be  sepa- 
rated. If  the  inverse  transformation  can  be  realized,  the  recovered 
signal  has  been  obtained  by  what  is  referred  to  as  generalized  linear 
filtering. 

The  theory  has  been  implemented  with  some  success  in  what  may  be 
called  multiplicative  systems.  If  x,  the  signal  to  be  processed, 

n  fi 

consists  of  components  x.  and  x^  combined  as  in:  x  =  x>  '  ^o   *    then  the 
transformation  employed  might  be  the  logarithm  function  (real  or  complex) . 
What  then  results  is, 

log  X  =  a -log  x.  +  P'log  Xg 
so  that  the  modified  components,  log  x^  and  log  x  ,  have  been  combined 
in  the  ordinary  linear  sense.  Whether  or  not  anything  has  been  gained, 
as  well  as  the  form  of  the  processing  to  follow,  depends  upon  certain 
properies  of  x.  and  x  ,  Another  class  of  systems  to  which  the  theory 
has  been  applied  may  be  called  convolutional. 

If  X,  the  signal  to  be  processed,  consists  of  components  x.  and 
Xp  which  are  combined  by  a  convolution,  denoted  x  =  x.  *  x  ,  the 
generalized  linear  filtering  technique  may  be  effective.  The  trans- 
formation might  consist  of,  first,  the  evaluation  of  the  Fourier 
Transform  since  it  maps  a  convolution  to  a  multiplication: 
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f[x]  =  F  [xj  -F  [x^] 


As  before,  this  can  now  be  followed  by  a  logarithm  function  evaluation. 
Here,  the  logarithm  must,  in  general,  be  complex.  What  results  is, 


log  F 


xJ  =  log 


^1 


+  log  F 


Again,  whether  or  not  anything  has  been  gained  depends  on  the 
properties  of  x.  and  x^,  Schafer    has  examined  a  specific  class 
of  signals  in  which  the  x>  and  x  exhibit  favorable  properties.  His 
presentation  is  closely  followed  herein.   The  signals  may  be  described 
as  echoed,  i.e.,  of  the  form 


(1.1) 


x(t)   =    La  s(t-T  ) 
1=0        ^ 


Such  a  signal  consists  of  a  basic  waveform  with  scaled,  delayed  versions 
of  the  same  waveform  superimposed.  This  is  a  reasonably  approximate 
description  of  multipath  signals  which  can  arise  in  many  applications. 
The  components  x>  and  x-  are  identified  by  re-writing  equation  (l.l) 
as  a  convolution  of  the  basic  waveform  and  an  impulse  train, 


(1.2) 


x(t)  =  x^(t)  *  X2(t) 
x^(t)  =  s(t) 


X  (t)  =    La  6(t-T.) 
'^  i=0 

where  6(')  represents  the  Dirac  delta  function 


The  important  aspect  of  the  identification  made  in  equation  (1.2) 


is  that 


log  F 


log  F  [x^ 
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is  a  slowly  varying  function  of  frequency  while 


is,  in  many  cases,  a  rapidly  Vcorying  function.  Consequently, 


the  inverse  Fourier  Transforms  of  log  F 


L^lJ 


and  log  F 


L^2 


are,  in 


many  cases,  somewhat  disjoint.  A  more  important  property  is  that 


F~^  log  F  [x^' 
„-l  ' 


is  also  an  impulse  train.  Hence,  by  use  of  the 
log  F  I  X  I  ,  the  convolved  components  are  transformed 


F-1 

log 

F  [x" 

operation  F 

to  added  components  which  exhibit  properties  that  may  make  separation 

possible. 

A  similar  operation  was  described  by  Bogert,  Healy  and  Tukey 
in  which  the  power  spectrum  of  the  log  of  the  power  spectrum  of  the 
convolved  signals  was  used.   The  result  was  called  the  cepstrum  of 
the  signal,  denoted  x, 


A 

X  = 


Nany  properties  of  both  operations  are  similar.   In  particular, 
the  impulse  train  is  mapped  to  another  impulse  train  in  both  methods. 
A  major  difference  is  that,  with  the  cepstrum,  the  loss  of  the  phase 
information  precludes  any  inverse  transformation.  Although  filtering 
and  recovery  ai-e  not  possible,  some  parameters  of  the  signal  are 
measurable  in  the  cepstrum  so  that  the  technique  finds  many  appli- 
cations. Due  to  the  similarities  involved,  the  transformation  which 
retains  the  phase  information  has  come  to  be  known  as  the  complex 
cepstrum  while  the  other  is  sometimes  referred  to  as  the  power 
cepstrum  to  maintain  the  distinction. 

On  the  subject  of  terminology;  it  becomes  awkward  to  speak  of 
"rapidly  varying  functions  of  frequency,"  as  above,  when  the 
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temptation  is  to  say  "high  frequency  components  -  in  the  frequency 
domain."  In  this  regard,  some  phrases  have  been  developed  in  the 
literature  of  this  subject  and  will  be  used  to  avoid  confusion.   In 
this  thesis,  the  term  power  cepstrum  will  be  used  to  refer  to  the 
transformation  of  Bogert  et  al.  The  transformation  developed  by 
Schafer  will  be  referred  to  as  the  complex  cepstrum,  or  as  simply 
the  cepstrum.  The  cepstinm  will  be  considered  to  be  a  function  of 
quefrency.  Thus,  if  the  complex  logaxithra  of  the  Fourier  Transform 
of  a  signal  has  a  component  which  is  a  rapidly  varying  function  of 
frequency,  it  has  a  high  quefrency  component.  Filtering  operations 
on  the  cepstrum  will  be  referred  to  as  short  pass  (analogous  to 
low  pass  but  in  the  quefrency  domain) ,  long  pass  (analogous  to  high 
pass),  and  comb  (the  usual  -  but  in  the  quefrency  domain)  operations. 

With  the  advances  in  technology  leading  to  high  speed  digital 
computation  and  fast  a/d  conversion,  and,  with  the  development  of 
the  Gooley-Tukey  algorithm,  realization  of  the  various  mathematical 
operations  called  for  by  the  technique  has  become  increasingly 
practical.  Since  the  implementation  is  via  digital  computer,  the 
theory  has  been  developed  in  terms  of  sampled  data  versions  of  the 
convolved  signals.  An  in-depth  presentation  of  the  theory  is 
contained  in  Schafer   and  only  the  relevant  portions  will  be  related 
herein.  A  brief  overview  of  transform  theory  for  sampled  data  signals 
is  called  for  before  the  specifics  are  presented, 

A  continuous  -  time  signal,  f(t),  may  be  sampled  every  DT  time 
units  to  produce  a  sequence,  f(n)  =  f(n-DT).  Such  a  sequence  has  a 
Z  -  Transform  defined  as. 


-2&. 

-n 


F(z)  =    E   f(n)  z 


The  inverse  operation  is  accomplished  via: 

f(n)  =  ■—    §^   F(z)  z"-l  dz 

where  G  is  some  closed  contour  in  the 
region  of  convergence  of  F(z) 

By  this  operation,  f(n),  a  function  of  a  discrete  variable  is 
transformed  into  a  function  of  a  continuous  complex  variable.   The 
Fourier  Transform  of  a  discrete  sequence  is  defined  as  the  Z-Transforra 
evaluated  around  the  unit  circle  in  the  z  -  plane, 

F(eJ")  =    E   f(n)  e'^"" 
.  n=-«> 

The  inverse  operation  is, 

f (n)  =  ^  1     F(eJ")  eJ^  dw 

-TT 

Both  the  Z-Transform  and  the  Fourier  Transform  will  be  referred 

to  in  the  development  of  the  report  since  they  allow  easily  accomplished 

evaluations  of  results.  Neither,  however,  can  be  directly  realized 

by  a  digital  computer  since  they  are  functions  of  continuous  variables. 

What  is  utilized  is  the  so-called  Discrete  Fourier  Transform  (DFT) 

which  is  the  Z-Transform  evaluated  at  equally  spaced  points  axound 

the  unit  circle  in  the  z  -  plane, 

N-1        .i2-r^n  ' 
F(k)  =    E  f(n)  e"   N 
n=0 
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and  the  associated  inverse  is 

.       N-1  .i2nkn 

f(n)      =     -~-      L      F(k)      e       N 
"       k=0 

All  of  these  transforms  may  suffer  from  shortcomings  referred  to 
as  aliasing  and  leakage.  Additionally,  the  DFT  is  accompanied  by 
the  so-called  picket  fence  effect.   These,  effects  and  means  of 
counter-acting  them  are  developed  by  Gcoley  et  al  ,  and  will  be 
discussed  in  this  report  only  when  they  present  problems.   The  DFT 
is  evaluated  via  the  Cooley-Tukey  algorithm,  or  the  Fast  Fourier 
Transform  (FFT), 

In  the  case  of  interest,  the  Z-Transform  of  the  sampled  composite 
signal  is, 

X(z)  ^    X^(z)-X^(z) 

The  logarithm  is  taken  to  obtain, 

X(z)  =  log  X(z)  =  log  X^(z)  +  log  X^Cz) 

The   complex  cepstrum,  "xCn),  is  to  be, 

:(n)  =  Z-^  [x(z)]   =  2^  #^  X(z)  z""^  dz 


x{ 


Here  the  complex  logarithm  presents  a  problem  since  it  branches 

A 

at  values  of  z  for  which  ARG  X(z)   =  +tt.  When  this  happens,  X(z) 
has  a  discontinuous  imaginary  part  and  is  therefore  not  analytic. 
Since  the  cepstrum  will  be  computed  via  the  DFT,  it  is  required  that 
ARG  X(z)  be  continuous  for  values  of  z  along  the  unit  circle.  The 
problem  is  overcome  via  a  so-called  phase  unwrapping  technique  which 
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identlfies  the  proper  branch  of  the  arctangent  function  defining 
the  phase.  This  technique  is  described  in  Appendix  B. 

More  specifics  of  the  computation  will  be  presented  in  subsequent 
chapters.  The  important  aspects  of  the  technique  can  be  siimmarizedj 

1.  Convolutions  are  mapped  to  additions, 

2.  The  mapping  is  continuous;  the  inverse  mapping  exists 
and  is  continuous  -  hence  the  term  homomorphic 
deconvolution, 

3.  Deconvolution  is  accomplished  via  a  subtraction  (if 
it  is  known  what  to  subtract) , 

4.  The  components  of  the  original  signal,  in  some  cases, 
become  somewhat  distinctive  in  the  complex  cepstrura 
so  that  the  may  be  identified  and  separated. 

The  technique  has  found  many  applications  since  it  can  deconvolve 
signals  when  the  components  are  unknown.   It  is  merely  necessary  that  the 
signal  be  reasonably  echoed  in  nature  with  reasonably  spaced  echo 
arrival  times.  The  basic  waveform,  either  as  it  is  or  with  minor 
processing,  will  have  a  low  quefrency  nature  and  decrease  in 
magnitude  at  least  as  fast  as  l/n  in  the  cepstriim.   The  impulse  train 
can  be  made  to  have  its  lowest  quefrency  component  displaced  from  the 
zero  quefrency  line  by  a  value  corresponding  to  the  first  echo  arrival 
time.   In  cases  of  large  echo  times,  the  echoes  can  be  removed  by 
a  simple  short  pass  or  more  involved  comb  filter.  Conversely,  the 
basic  waveform  can  be  removed  and  the  impulse  train  recovered  by 
use.  of  a  long  pass  filter.  This  approach  is  applied  in  cases  in 
which  the  multipath  response,  x„,  is  not  strictly  impulsive  and  it 
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is  desired  to  see  precisely  what  it  is.   In  such  a  fashion,  properties 
of  the  paths  involved  can  be  measured. 

Some  applications  of  the  technique  have  been  to  speech  proces- 

3 
sing  by  Schafer,  Oppenheim  and  Stockham  ,  to  seismology,  where  it  is 

7 
desired  to  recover  and  measure  x  ,  by  Ulrich  ,  and  to  processing  of 

measured  brain  waves  evoked  by  visual  stimuli  by  Senmoto  and  Childers  . 

In  this  last  application,  the  data  is  available  via  an  electro- 
encephalogram. There  is  some  degree  of  control  in  that  the  echoed 
signal  measured  represents  the  composite  response  of  the  brain  to  a 
pulse-like  visual  stimulus.   The  stimulus  can  be  re-applied,  at 
fixed  intervals,  over  a  pulse  repitition  interval.  A  somewhat 
involved  technique  is  used  to  average  the  responses  obtained  over 
a  repitition  interval  so  that  improved  signal  to  noise  ratio  is 
obtained  before  computing  the  cepstrum. 

The  similarities  between  this  application  and  the  Loran-G 
case  lead  to  speculation  on  signal  to  noise  ratio  improvement 
which  will  be  discussed  in  Chapter  6,  Considerations  such  as  the 
effects  of  interference,  atmospheric  noise,  and  signal  stability 
must  be  taken  into  account.  Before  these  are  considered  however, 
more  specifics  of  the  problem  must  be  developed. 


-32- 
Ghapter  2 
PROBLEM  FORMULATION 

In  order  to  apply  the  homomorphic  deconvolution  technique,  a 
model  of  the  Loran-G  signal  as  described  in  Chapter  1  must  be 
developed.  From  this  model  the  signal  processing  tasks  can  be  identi- 
fied so  that  an  appropriate  solution  scheme  results.   To  this  end, 
the  waveform  to  be  processed  will  be  considez-ed  to  be  obtained 
from  one  pulse  duration  of  the  received  signal,  represented  as 
follows: 


k  f- 

(2.1)        F(t)  =  n(t)  +    S  A,  s.(t)  sin  w^(t-T,)  -9, 

i=0  ^  ^       L  c    1    1 


In  equation  (2.1),  the  zero  subscripts  identify  parameters  to 
be  associated  with  the  received  groundwave.  The  subscripts  1  through 
k  identify  parameters  of  the  first  through  k   skywaves  respectively. 
s. (')  denotes  the  normalized  pulse  envelope  of  the  i   component 
of  the  composite  signal.  The  T. 's  and  A. 's  refer  to  the  envelope 
arrival  times  and  magnitudes,  respectively,   9.  represents  the 
phase  shift  between  the  envelope  and  the  carrier  of  the  i 
component. 

With  a  knoi-m  carrier  frequency,  an  equivalent  method  of  speci- 
fying the  signal  is  by  use  of  the  complex  envelope  representation 
cis  in  equation  (2.2). 


(2.2)       *(^^  ^  *^^)  ^    2  .   g^p 

i=0  ^ 


w  T.  +9, 
c  i    i 


.s.(t) 
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Figure  2-1  Generation  of  Complex  Envelope 

A  common  method  of  describing  how  the  complex  envelope  is 
obtained  is  shown  in  figure  2-1.  An  analysis  of  the  signal  flow 
in  the  figure  shows  that, 


f  (t)  =  n(t)  cos  w  t  -    E  A,  s  (t)  sin  w  T.  +  0. 

C  C         ._^   11  LCI     1 


and, 


f  (t)  =  n(t)  sin  w  t  +    E  A.  s. (t)  cos 

•  s  C         ,^11 

1=0 


w  T.  +  e, 

C  1     1 


and 


If  the  following  definitions  are  made: 


n(t)  =  n(t)  sin  w  t  -  j  n(t)  cos  w  t 


f(t)   =  fg(t)   -  j  f Jt) 


it  follows  that, 


(2.3)    f(t)  =  n(t)  +    E  A  exp 

i=0  ^ 


w  T.  +  e^ 
c  i    i 


Si(t) 
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which  is  identical  to  equation  (2.2).  Clearly,  the  complex  envelope 
representation  for  continuous  signals  is  no  more  than  a  mathematical 
convenience.   In  a  sampled  data  form  however,  it  can  be  actually 
implemented  and  provides  a  compact  means  of  handling  the  quadrature 
components  of  the  baseband  signal.  The  realization  method  is  as 
depicted  in  figiire  2-2. 
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f  (n-bT)  =  f^(n) 


R' 


90' 


Lovf  Pass 
Filter 


Sampler 
(DT  sec) 


f  (n.DTl_  fjCn) 


Analog  I  Digital 


Figure  2-2  Realization  of  Sampled  Complex  Envelope 


The  component  f  (t)  is  sampled  and  stored  in  the  real  part  of 
the  sample  word  while  f  (t)  is  sampled,  complemented,  and  stored  in 
the  imaginary  part.  What  results  is  a  sequence  of  the  following 
form» 


f(n)  =  n(n)  +    S  A  exp 

i=0  ^ 


J  \\\   +  G. 


s.  (n) 


For  the  remainder  of  the  presentation  in  this  chapter,  the  signal 
will  be  considered  to  be  noise  free.  The  theory  of  noisy  signals 
is  developed  in  Chapter  5«  Consequently,  the  signal  of  interest  is, 
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f(t)  = 


i;   A.  exp 
i=0  ^ 


Vi-^^i 


s^Ct) 


For  notational  convenience,  and  without  loss  of  generality,  the 
magnitudes  will  be  normalized  with  respect  to  A^  and  the  phase 
notation  will  be  simplified: 


A'      =     — 
i  A, 


t       1—U|     1|     ••«|K 


'^i  oil 


B.      =     A. 
1  1 


.'    e'^^i 


Hence, 

(2.4)  f(t)     = 


2     B     s   (t) 
i=0     ^     ^ 


In  Chapter  4,  distorted  signals  will  be  considered.  Until  then, 
it  will  be  assumed  that  s. (t)  =  s(t-T. )•  In  this  case,  equation 
(2.4)  becomes 


f(t)  =    S  B  s(t-T  ) 
i=0  ^      ^ 


This  is  "Uie  form  of  signals  for  which  the  horaomorphic  deconvolu- 
tion  technique  has  been  developed.  Additionally,  it  is  possible  to 
identify  the  tasks  to  be  performed  at  this  point.  As  presented  in 
Chapter  1,  it  is  desirable  to  measure  0^,      the  groundwave  parameter 
from  which  time  measurements  are  made,  and  T^,  to  measure  the  difference 
in  propagation  speeds  of  the  envelope  and  carrier.  Additionally, 
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measurements  of  T^  and  B^   would  be  helpful.  Finally,  measurements 
of  all  ^. 's,  T, 's,  and  B,  's  would  be  desirable  since  a  total 
power  receiver  might  then  be  feasible. 

With  the  problem  thus  defined,  solution  techniques  can  be 
considered.  To  begin,  consider  the  following  to  be  a  reasonable 
approximation  to  the  received  Loran-C  pulse  envelope: 


.(t) 


,2  -at    /,\ 
t  e    u  .(t; 


a  =  ("ZP-;  X  10   sec 


where  u  ^ (t)  is  the  Heaviside  unit  step  function 


Initially,  make  the  following  assumptions: 


K  =  ° 


T.   =  n  •  DT 
1      i 


n.   =  an  integer 


DT  =  the  sampling  time 


The  sampled  waveform  in  this  case  is; 


k  k 

f(n)  =    Z  B.  s(n-n.)    =  s(n)  (g)    E  B.  6(n-n.) 


i=0 


i=0 


where  (g)  denotes  a  discrete  convolution  and  6(n-n. )  is  a  unit  sample: 


6(n-n^)  = 


n  =  n. 


n  5^  n. 
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Therefore, 


where 


f(n)  =  s(n)  ®  p(n) 


k 
p(n)   =    E  B  6(n-n  ) 

In  the  case  of  k  =  1,  I.e.,  in  the  case  of  a  signal  consisting  of 
the  groundwave  and  one  skywave, 

p(n)      =  <5{ti)+  B^  6(n-n^) 


If    |B^ 


<   1,    the  complex  cepstrum  of  p(n)    is      (see  Appendix  A  for   the 
details  of  the  computation) : 

$(n)   =    Z  (-1)'"^^  -^  i(n-mn^) 

m=l 
whereas  the  cepstrum  of  s(n)  is, 

(2.5)  ^(n)  =  -<x6(n)  +    E  ^  "^  '^"^^"'    e"°^'"6(n-m) 

m=l 


Since  the  complex  cepstrum  operation  has  the  property  of  mapping 
convolutions  in  the  time  domain  to  additions  in  the  quefrency  domain, 
the  composite  cepstrum  consists  of  s(n)  +  p(n).   In  this  case,  a 
simple  estimator  can  be  made  to  scan  the  cepstrum  sequence  and 
determine  n>  from  the  location  of  the  components  of  the  impulse  train. 
A  simple  comb  filter  could  then  be  employed  to  remove  all  traces  of 
the  echo  from  the  cepstrum  -  only  slightly  effecting  the  component 
due  to  s(n).  Before  filtering  p(n) ,  an  estimate  of  |B||  and  0^ 
could  be  made  from  the  components  of  the  impulse  train.   In  such  a 
manner,  estimates  of  T^ ,  |B.  I,  and  ^i  ^-^^  obtained  and  all  skywave 
contamination  is  removed. 
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When  T.  is  not  an  integer  number  of  sampling  times  greater  than 
T_,  ^(n)  no  longer  has  such  a  simple  form.  The  components  of  the 

impulse  train  become  terms  which  spread  throughout  the  cepstrum 

sequence  as  in  equation  (2.6)  (See  Appendix  A), 


(2.6) 


oo 


r,in 


P(n) 


™       ,1   B.    sin 
2  (.1)^-^1  _1 


Tr[n-m(n^+A^); 


n-m(n^+A^)J 


where     n.   =  an  integer 
and      T^     =  DT'(n^+A^) 
-  0.5  <  Aj^  <  0.5 
In  this  case  the  filtering  becomes  less  straightforward,  A  simple 
comb  filter  can  not  be  effective.  In  the  distortion  free,  noiseless 
case  however,  the  cepstrum  consists  only  of  s(n),  which  is  known, 
and  p(n),  which  is  known  to  within  three  parameters:  B.  ,  n.,  andAi» 
so  that  an  estimator  can  be  developed  to  identify  p(n).  When  this 
identification  is  made,  the  deconvolution  can  be  accomplished  via 
a'  subtraction.  The  form  of  the  estimator  is  developed  in  Chapter  3 
and  the  effects  of  signal  distortion  and  noise  on  its  performance 
are  presented  in  Chapters  4  and  5« 

The  problem  is  further  complicated  when  T^   is  not  an  integer 
multiple  of  the  sampling  timoo   In  this  case,  the  s(n)  of  equation 
(2,5)  has  another  term  added  to  it,   (Appendix  a) 

vn+1   Ao 


(2.7) 


s'(n)  = 


A 


s(n)  +  (-1)^ 


^(n) 


n 


;   n  ^  0 


n  =  0 
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The  presence  of  this  new  term  makes  the  estimation  procedure  described 
above  suspect.  Now,  there  are  two  unknown  components  in  the 
computed  cepstrum.  Since  the  forms  of  the  components  are  known 
however,  a  modified  estimator  can  still  be  effective. 

A  summary  of  the  problem  as  formulated  thus  far  can  be  made 
via  figure  2-3 •  The  received  signal  is  processed  to  produce 
quadrature  baseband  components.  These  are  sampled  and  digitally 
combined  to  produce  the  sampled  complex  envelope.   The  complex 
cepstrum  is  computed.   From  the  cepstrum,  signal  parameters  are  estimated. 
One  of  the  uses  of  the  estimated  parameters  is  shown:  the  groundwave 
phase  estimate  provides  feedback  to  the  local  oscillator  to  permit 
phase  tracking.  An  important  aspect  of  this  particular  use  i s  that 
it  provides  a  method  of  comparing  the  performance  of  the  homomorphic 
filter  with  other  receiver  techniques. 
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Ghapter  3 
DEJCONVOLUTION  OF  NOISELESS  SIGNALS 

As  developed  in  Chapter  2,  it  is  desired  to  estimate  vraveform 
parameters  from  the  cepstrum  of  the  received  signal  so  that  the 
separation,  or  deconvolution,  can  be  achieved.  In  the  distortion 
free,  noiseless  case,  the  technique  is  as  presented  in  this  Chapter. 

It  is  presumed  that,  with  one  skywave,  the  received  signal  will 
have  a  cepstrujn  of  the  following  form: 


f(n) 


(3.1) 


s'(n)  +  9(n) 


=  s(n)   +  g(n)   +  p(n) 


where,  from  equations  (2.5) »  (2.6),  and  (2.7), 


m+1 


(a)    a(n)  =  -a  6(n)  +  E  ^-^-^^ e"^™  6  (n-m) 

Iin 


(3.2) 


(b)    g(n)  = 


0  ;   n  =  0 


/  ^Nn+1  _A0    ;   n  7^  0 


n 


(c)    ^(n)  = 


CO 

z  (-1) 

m=l 


m+1   ^1 


m 


sm 


ra 


Tr[n-m(n^-+A^)j 
TT[n-m(n^+/!^j^)J 


In  the  estimation  procedure  it  is  desired  to  use  the  knowledge  of 
the  form  of  the  above  sequences  in  processing  the  received  signal. 
Consequently,  it  will  be  necessary  in  the  development  to  refer  to  the 
above  sequences  and  corresponding  components  of  the  received  signal 
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simultaneously.  For  clarity  of  presentation  therefore,  the  notations 
s(n)i  6(n),  p(n),  and  f(n)  will  be  used  to  refer  to  the  predicted,  or 
ideal,  sequences  as  defined  above  whereas  the  received  signal  will 
be  described  as, 

r(n)  =  s  (n)  ®  g  (n)  ®  p  (n) 

(3.3)                 a  a        a 
or 

r(n)  =  s^(n)  +  g^(n)  +  p^(n) 

The  subscript  a  signifies  actual  components  of  the  signal  as 
received.  In  the  distortion  free  case,  the  only  difference  between 
the  components  of  f(n)  in  equation  (3«2)  and  r(n)  in  equation  (3.3) 
is  due  to  computation  error, 

A  third  set  of  sequences  will  also  be  referred  to.  Estimates  of 
signal  parameters  Aq»  Ah »  n. ,  and  B^  will  be  made  from  the  received 
signal.  From  these,  estimated  sequences  will  be  generated.  The 
estimated  sequences  will  have  the  form  of  those  in  equation  (3.2) , 
but  with  parameters  as  estimated  from  the  sequences  of  equation  (3.3) o 
Thus,  the  estimated  sequence  of  g(n)  is: 


g(n)  =  g(n) 


Aq  =Ao 


(   iNn+l  Ao         /  _ 
(-1)    —^       !   n  ^  0 

0         I   n  =  0 


where  Aq  is  the  estimate  of  Af^  obtained  from  r(n). 

3.1  Sequence  of  Estimation 

The  received  signal  consists  of  the  terms  s  (n) ,  g  (n),  and  p  (n) 
which  are  combined  to  produce  the  cepstrum  sequence  shown  in  figure 
3-1. 


-43- 


'  ,  '  I 


-1 


-2 


H 1 1 1 1 i 1- 


-5 


H 1 ^-+- 


Figure  3-1  Gepstrum  of  Simulated  Received  Signal 


-i4h. 
It  can  be  seen  that  there  is  insufficient  seperation  between  the 
sequences  to  allow  for  precise  estimation  of  the  various  parameters 
directly.  Consequently,  the  estimation  is  accomplished  in  a  number 
of  steps;  systematically  estimating  the  parameters  one  at  a  time. 
The  procedure  used  herein  is  to  make  the  estimates  in  the  following 
order. 

1.  ^(n)  -  an  estimate  of  the  received  signal  basic 

waveform. 

2.  ^(n)  -  an  estimate  of  the  small  linear  phase  term 

not  removed  before  computing  the  cepstnira 

3.  ^(n)  - 

(a)  n^   -  a  coarse  measurement  of  the 

skywave  arrival  time 

(b)  A-i   ~  3-  precise  measurement  of  the 

skywave  arrival  time 

(c)  B>   -  skywave  magnitude  and  phase 

The  procedure  begins  by  utilizing  the  fact  that  a  good  a  priori 

estimate  of  s  (n)  is  available,  to  wit:  s(n)  of  equation  (3.2). 
a 

Therefore  the  first  estimate  is  s(n)  =  s(n).   This  estimated 
sequence  can  be  subtracted  from  the  received  cepstrum  to  generate  a 
new  sequence  denoted  x(n;, 


x(n)  =  r(n)  -  s(n)  =  g^(n)  +  p^(n)  +  € ^(n) 


where 

6  (n)  =  s  (n)  -  s 


s        a 


(n) 


In  this  chapter  it  can  be  assumed  that  this  error  sequence,  f  (n;, 

is  negligibly  small.   The  procedure  as  thus  far  discussed  was  accomplished 

on  a  simulated  received  signal  and  the  resulting  x(n)  is  shovm  in 
figure  3-2. 


Figure  3-2  Simulated  x(n)  Sequence 


It  can  be  seen  that  there  is  improvement  over  the  situation  shovm  in 

figure  3-1 »  but  that  the  components  of  g  (n)  and  p  (n)  still  interfere 

a        a 

so  that  the  desired  estimates  are  still  not  readily  obtained.  More 
precisely,  it  can  be  seen  that  for  quefrency  values  greater  than 


(n^^+Ai),  the  various 


sin  X 


terms  of  ^  (n)  interfere  with  one  another 


to  some  extent.  For  positive  quefrencies  less  than  (n.+A-),  there  is 


-46- 

interference  between  g  (n)  and  the  first  term  of  p  (n). 

What  makes  the  estimation  possible  is  the  fact  that  for  small 
negative  quefrencies,  p  (n)  is  reduced  enough  so  that  it  only  slightly 

Si 

interferes  with  g  (n).   The  next  step  in  the  procedure  therefore 
will  be  to  estimate  Af^  from  the  negative  quefrency  portion  of  the 
cepstrum. 


3»2  Estimation  of  g  (n) 


Since  |(n)  =  (-l)""'^  _Ao  .  „  <  q, 


it  follows  that 


Ao  -    i(-i) 


-^°-  =  l(-2) 


2      "^^a^ 
—  -  g^(-3) 


where  A^v  will  be  the  estimate  of  An*  As  described  above,  it  is  true 

that,  - 

r(-l)  =  ^^(-1)  +  ^^(-1)  +  ^^(-l)  «   g^(-l) 

r(-2)  =  g^(-2)  +  ^^(-2)  +  e^(-2)  «  g^(-2) 


since  the  p  (n)  sequence  is  small  for  this  range  of  quefrencies.  The 

3. 

estimate  of  A^  ca-n  therefore  be  obtained  from, 


-/I?. 


(3.^) 


1      1      1 


^, 


With  Aa  so  estimated,  the  sequence  g(n)  can  be  generated  and 
subtracted  from  x(n)   to  produce  the  new  sequence  y(n): 


y(n) 


=  x(n)  -  g(n) 


where 


Kn)  = 


(-1) 


n+1   Ao 
n 


n  =  0 


n  f^  0 


The  procedure  was  applied  to  the  x(n)  in  figure  3-2  to  produce 
the  y(n)  of  figure  3-3. 


A  ,   . 

Figure  3-3  Simulated  y{n)   Sequence 


-48- 


At  this  point  in  the  procedure,  since  T-  =  DT  (n^+ls^) ,    the 


estimate  of  T„  can  be  made: 


%     =  I^T  (ri^^/^) 


where  n.  is  deduced  in  the  phase  unwrapping  computation  as  described 
in  Appendix  G, 

The  reason  for  the  particular  estimator  used  in  equation  (3.^+) 
can  be  seen  from  figiire  J-k, 


-1012345678 


Figure  3-4     Comparison  of  g(n)  and  p(n) 
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Both  sequences  g  (n)  and  p  (n)  decrease  essentially  as  fast  as  l/n 

3.  3. 

for  negative  values  of  n.  Near  the  zero  quefrency  value,  the  g  (n) 
terra  clearly  dominates  the  composite  signal.  As  n  becomes  more 
negative  however,  the  distinction  between  the  two  components  is  not 
as  great.  Consequently,  only  a  few  samples  should  be  used  to  make 
the  estimate  and  even  among  these  few,  the  samples  obtained  from  the 
more  negative  quefrencies  should  be  less  heavily  weighted.  This  is 
precisely  what  is  accomplished  in  equation  (3.^), 

The  next  step  in  the  procedure  is  to  estimate  the  parameters 
of  p(n),  i,e,,  All  n.  ,  and  B^ ,  from  the  sequence  y(n).  This  is 
accomplished  in  three  steps,  beginning  with  the  estimation  of  n^. 

3,3  Estimation  of  p  (n) 

The  values  of  y(n)  are  examined  for  quefrencies  corresponding 
to  times  between  about  30  microseconds  and  60  microseconds.  The 

largest  value  of  y(n)  in  this  range  of  quefrencies  will  be  assumed  to 

sin  X 
correspond  to  a  sample  in  the  main  lobe  of  the  first term  in 

p  (n).  Consequently,  the  value  of  n  at  v;hich  this  maximiim  occurs 

will  be  the  estimate  of  n.  ,  The  situation  would  typically  be  as 

depicted  in  figure  3-5»  Since  the  values  of  y(n)  will,  in  general, 

be  complex,  the  plot  in  the  figure  corresponds  to  either  the  real 

or  iminaginary  part  of  y(n). 
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Figure  3-5  Expected  Form  of  Pa(n) 

Once  n^  is  determined,  B^  and  A-t   can  "be  estimated  as  follows, 
Consider  the  values  of  y(n)  in  the  vicinity  of  n  =  n.   to  be  due 
only  to  the  first 


function  in  the  expansion  of  p(n).  As  an 


example  of  the  validity  of  this  assumption,  consider  the  values  of 
the  first  and  second  terms  in  the  sum  of  equation  (3.2b)  for  n  =  n^ : 


p(ni) 


sin  Trfn^   -   (n^+A^)]  B^       sin  Trfn^  -  2(n^+A^) 

iT[n^  -   (V-^l^j  "     "^         TTfn^  -  2(n^+Ai)J 


+  .    .    . 


B>    sin  ttA. 

a; 


B^^sin  TT  [(n^  -  2A^)J 

2'rT-(n^-2Aj^)  "*"  *    * 


The  ratio  of  the  magnitude  of  the  second  term  to  the  magnitude  of  the 
first  term  will  be, 
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|Al|    |Bi|      l^i"  ^Ail 

2    |n^  -2A^|    jB^I    isin-rrA^I 


Ai 


n^  -  2A^ 


sin  2rTA^ 


2  sin  ttAj^ 


Exponential  weighting  ,   as  described  in  Appendix  A,   ensures  that 
less  than  1.     Additionally,    the  term 


is 


sin  2ttAi 


2  sin  TiZ\. 


is  less  than  1  so  that  the  ratio  is  less  than 


(3.5) 


n^  -  2A^|  2  n^  -  ^A^ 


for     n.  >    1 
since   |AJ   <  f 


Hence,    for  reaLSonably  large  values  of  n,  ,    the  approximation 
is  a  good  one     and  allows  the  following: 

B^  sin  irfn^  -  2  -   (n^+A^)]  B^  sin  Tr(2+A^) 


y(n^-2) 


TTJ^n^  -  2  -   (n^+A^)J 


•n(2+A^) 


B.    sin  Tifn.    -1  -   (n.+A.  )'!                 B.    sin  tt(1+A,) 
/^f         jk\      _        i L   -^ 1      1  J  _  1 i. 

y^"r^^     -  ^|n^  -  1  -   (n^Hii^)j  -  -nd+A^ 


Bj  sinn[n^  -   (n^+Apl 
(3-6)        y(ni)  =    TTfn^  -   (n^-fA,)] 


B.    sin  ttA. 


"ST 


■^                        B^  sin  TT^n^  +  1  -   (n^+A^)!               B^  sin  Tr(l+A^) 
^^"l""^^      ^    nfn^  +  1  -   (n^+A^Vj =  rrd+A^) 


Bj^  sin  TT[n^  +  2  -   (n^+A^)!  B^  sin  tt  (2-+A-||) 


^^'^l''^^      =     TTf^a^  ;  2  -   (n^+Ai)j 


tt(2+A^) 


The  solution  can  proceed  in  two  steps:   estimate  A< ,    then     B^.     To 
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accomplish  this,  define  the  following  variables! 
^  y(n^-2)        Al 


^1 


(3.7) 


^2  =  -J(^ 


d 


y(n^+i) 


^k 


3  "   y(n^) 


2 

+  Ai 

-Al 

1 

■^Ai 

Al 

1 

-Al 

-Al 

2 

-Al 

It  can  be  seen  that  fovir  estimates  of  A^  can  be  obtained  from 
the  d. 's.  The  estimates  obtained  from  d  and  d^  however,  are  made  from 
values  of  y(n)  which  are  closer  to,  or  in,  the  main  lobe  of  the  first 

—  term.  They  are  therefore  less  sensitive  to  the  effects  of 

interference  from  other  terms  and  the  residuals  of  the  "g  (n) 

estimate  than  the  estimates  obtained  from  d.  and  d^.   It  is,  however, 
desirable  to  base  the  estimate  on  as  many  measurements  as  possible 
so  that  a  compromise  would  be  to  maice  a  final  estimate  based  on  a 
weighted  sura  of  the  four  individual  estimates.  The  procedure  is: 

1,  Compute  the  numbers  d^  ,  d  ,  d„,  and  d^^  as  in  equation 

(3.7) 

2,  Make  four  individual  estimates  as  follows: 
2d, 


1 


11      ^  -  '^i 


A 


12     ^  -^  ^2 
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^13 

= 

S 

l  +  d3 

0! 

^14 

= 

-"4 
1-d^ 

where  A^T  denotes  the  estimate  of  A^  based  on  the 

measurements  which  generated  d,  . 
3»  Weight  these  estimates  by  the  corresponding  magnitudes 
of  the  measiorements  generating  them  and  perform  an 
appropriate  scaling  to  ajrrive  at  the  final  estimate: 


2 


^        k=l 


\\^ 


Ik 


1  ~    4 

k=i  '  ^' 

When  this  estimate  of  Ai  is  obtained,  equation  (3.6)  can  be 
re-applied  to  estimate  B> .  To  do  this,  make  the  following  individual 
estimates, 

'A 

^^     sin  -ftCI  +A^) 

13  sin  'nAj_ 

tt(1  -  Ap    ^ 

14  =  sinnd  -A^)  ^^"l""^^ 


1.  =   -4^   Kn,) 


^15  ==  sin  ^(2  -  A^)  ^^"l'-^^ 
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These  are  combined  to  produce  the  final  estimates 

5 


Hn^   +  3  -  k)  I   B^j^ 


'1  ~    5 
2 

k-=l 


y(n^  +  3  -  k) 


with  the  estimates  of  n.  ,  A-i»  ^^^   B-i  »  "the  estimated  sequence 
p(n)  can  be  generated.  This  and  g(n)  can  be  subtracted  from  the 

A' 

received  signal  cepstrum,  r(n),  to  produce  s(n).  If  the  Inverse 
DFT  of  s(n)  is  computed,  what  would  result  would  be  an  estimate  of 
the  complex  log  of  the  groundwave.  From  the  zero  frequency  point 
in  this  sequence  can  be  measured  0^     -   the  phase  of  the  groundwave. 

If  there  are  more  than  one  skywaves  present  in  the  received 
waveform,  ^  (n)  will  have  additional  —-^ terms  centered  at 

3.  X 

quefrencies  of  (n  +Ap)  f  2(n  H-Ag).  .  .  .  i  (n^+Zu)  .  o  .,  etc.  and 
at  various  combinations  of  (n.-iA^)  and  (n^+Ap)  as  described  in 

Appendix  A,  Once  the  terms  due  soley  to  the  first  skyi-iave  are 

sin  X 
removed  however,  the  lowest  quefrency  — term  that  remains  is 

due  only  to  the  second  skywave  so  that  the  above  estimation  procedure 

can  be  re-applied  to  estimate  n^,  Ap»  and  B_.  This  can  be  done  as 

often  as  necessary.  The  entire  estimator  can  be  described  as  in 

figure  3-6. 
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3.^<-  Performance  of  Estimators 


Before  the  results  of  applications  of  this  technique  are  pre- 
sented, a  reason  for  the  ad  hoc  nature  of  the  estimators  can  be 
given.  Due  to  the  sequential  nature  of  the  estimation  scheme  and 
the  manner  in  which  the  received  signal  is  mapped  to  the  cepstrum, 
the  estimators  are  called  upon  to  act  in  the  presence  of  somev;hat 
structured  interference.  As  will  be  discussed,  the  interference 
structure  can  be  anticipated  so  that  the  estimators  can  be  designed 
to  provide  some  degree  of  protection  against  it.  To  show  this,  it 
should  not  be  necessary  to  work  with  the  complicated  estimators  as 
developed.   Instead,  a  similar,  but  much  simpler,  estimation 
situation  can  be  considered.  This  will  allow  a  clearer  presentation 
of  the  important  property  of  the  estimators. 

As  a  rough  model  of  the  sequence  to  which  the  A.  and  B.  estimators 
are  applied,  consider  the  sequence  z(n)  shown  in  figure  3-7. 


Figure  3-7  Simplified  Signal  to  be  Estimated 


;(n)  is  similar  to  p(n)  is  that  the  signs  of  the  z(n)  sequence 
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The  sequence  is  known  to  be  of  the  form  shown,  but  the  value  of  a 
is  unknown  and  to  be  estimated.   The  Important  manner  in  which 

z( 

values  have  a  partlculajr  code  to  them:  the  values  of  z(n)  are  positive 
for  odd  negative  values  of  n  and  even  positive  values  of  n.   This 

code  is  similar  to  that  encountered  in  p(n)  since,  in  p(n),  there 

sin  X 
are  two  samples  in  the  main  lobe  of  the  first  term,  causing 

a  similar  sign  reversal.  It  is  desired  to  consider  the  effects  of 

the  presence  of  structured  interference  on  the  estimation  of  a. 

To  do  this,  consider  first,  the  estimator  of  equation  (3*8)  which 

is  of  the  form  used  for  ^.    and  B> .   In  effect,  the  estimators  are 

matched  to,  or  Instep  with,  the  code  of  the  z(n)  sign  reversals. 


a 


1  =  -2\(-2) 


fo 


h     •  -a(-l) 


(3.8) 


a. 


A> 


-   -  ^.(1) 


Now  suppose  that  the  actual  sequence,   z  (n) ,   consists  of 

a 

the  anticipated  z(n)  as  in  figure  3-7i  plus  a  constant  interference 
sequence,  w(n)  =  b.  The  sign  code  of  w(n),  for  n  =  -2,  -1,  1,  2, 
16  +  +  +  +.  This  can  be  described  as  being  orthogonal  to  the  sign 
code  of  the  estimators,  which  is  +  -  +  -, 
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In  this  case,      z   (n)   =  z(n)  +  w(n)     =     z(n)   +  b,    so  that, 

3. 


^1 

c= 

-.(-i      ^ 

b) 

=     a  -   2b 

^2 

= 

(  a  +  b) 

=     a  +  b 

"3 

C3 

-   (  -a  +  b) 

=     a  -  b 

^4 

= 

2   (  -f-  +  b) 

=     a  +  2b 

Hence  the  individual  estimates  have  complementary  error.   It  should 
be  assumed  that  |a|  >2|b|  and  a>  0  in  making  the  following 
argument  about  the  combined  estimate. 


a  = 


-2)  I  ■  a^  +  |z^(-l)  |-  ^2  +  I  \(1)|-  a^  +  ra(2)|-  a^ 


Since     a  >  |2b|,   it  follows  that, 


\(-2) 

■■ 

'^          2 

+ 

b) 

(-f--b) 

\(-l) 

= 

(a  +  b) 

= 

(a  +b) 

\W     = 

(a-  b) 

= 

(a-b)      • 

>=,(2) 

= 

C      2     - 

b) 

= 

(-f-  +  b) 
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From  this,  it  also  follows  that, 

(I  _  b)  (a  -  2b)  +  (a  +  b)^  +  (a  -  b)^  +  (|  +  b)  (a  +  2b) 
(|-  b)  +(a  +  b)  +  (a  -  b)  +  (|  +  b)   • 

(3.9)        a  =  3^^ =  a  +  2-^ 


As  an  example,  suppose  that  b  =  a/5: 

a  =  a  +  TTP—  =  1.08  a 

The  estimate  for  this  technique  is  0.08a  whereas  if  only  z  (1), 
for  example, were  used,  the  error  would  have  been  b  =  0,2a,   If  b  =  ,1a, 
the  errors  are  0,02a  and  0,1a  respectively.  If  it  happens  that  the 
Interference  has  the  same  code  as  the  estimator,  it  can  be  shovm 
that  the  error  is  0.36a  for  b  =  ,2a  and  0,i7a  for  b  =  ,1a. 

If  H(n)  =  (-1)  b,  the  Individual  estimates  will  be, 

a^  =  -2(-  ~-  +  t>)   =  a  -  2b 


ap  =    (a  -  b)       =  a  -  b 


a^  =  -  (-a  -  b)      =  a  +  b 


aj^  =  2  (  -|-  +  b)    =  a  +  2b 


Once  again  the  errors  are  complementary  and,  it  can  be  shown,  the 
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total  estimate  will  be  as  in  equation  (3«9).  Figure  3-8  shows  the 
variation  of  the  estimate  error  as  a  function  of  b  for  an  interference 
sequence  of  constant  magnitude.  The  three  plots  are  of:  estimate 
from  z  (l)  alone,  four  step  estimate  -  interference  in  step  with 

Or 

estimators,  four  step  estimate  -  interference  orthogonal  to 
estimator. 


) 

/■    /' 

,5  a- 

/  y 

Interference    /                            y            / 
in  step  with  /                                          / 

^  pt  estimator                          /             / 
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/                       /              / 

^ 
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/                 y                     /orthogonal  to 
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/               y                      /  ^  pt  estimator 

/              "^                     / 

l/L"^                                                        1                                                                1 
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^      U 

' 

Interference  Level 

Figure  3-8  Estimator  Performance  Curve 
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In  the  actual  case,  the  estimators  are  a  bit  more  complicated 
than  those  of  equation  (3. 7)  and  the  symmetry  is  not  quite  complementary 
as  described.  The  example  does  however,  illustrate  the  improvement 

obtained  by  centering  the  measurements  around  the  main  lobe  of  the 

sin  X 

functions.  Two  cases  wherein  this  property  is  helpful  can 

be  presented. 

The  use  of  the  estimated  sequence  g(n)  will  generate  some 

small  residual: 

f  (n)  =  g^(n)  -  g(n) 

B  a. 

n      ^     ^  ^     ^     n 


This  is  one  source  of  interference  when  the  A^  and  B.  estimates  are  made. 
Although  the  sequence  is  not  of  constant  magnitude,  the  change  is 
only  fron  (^^„)/6  to  (fAf,)/lO  in  the  quefrency  range  of  interest 
for  n^  =  8,  This  is  a  reasonable  enough  approximation  to  a  constant 
magnitude  to  produce  some  improvement.  A  more  important  case,  since 
the  magnitudes  involved  will  typically  be  greater,  is  the  interference 

from  the  second  term.  As  developed  in  equation  (3«5)»  the 

magnitude  of  this  term  will  be  fairly  well  reduced  for  quefrencies 

near  n^ .  The  situation  can  be  examined  a  little  more  closely  to 

sin  X 
see  that,  if  the  largest  value  of  the  first  term  is  denoted  a, 

where  B^  sin  n^. 

a   =    -^— -r i  • 

TT  A^ 

then  the  interference  sequence,  denoted  w(n) ,  due  to  the  second 
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sin  X 


function  has  values, 


w(nj^+2) 
w(n^+l) 


B.   sin  ttA^ 


2iT(n^  +  A;!^  -  2) 

o 

B.   sin  rrZl, 


2i^(n^  +  A^  -  1) 


=   a  B 


=   a 


B. 


11  2(n^  +  Aj^  -  2) 
2(n^  +  A^  -  1) 


w(n^-2) 


B.   sinirA^ 
2Tr(n^  +Aj_  +  2) 


a  B 


l!  2(n^  +  A^  +  2) 


If  B^   -  0,8,  n^  =  8,  and  A^  -  'zt    these  values  range  from  about 
0,031a  to  about  0,019a.  Once  again,  the  symmetry  is  not  perfect, 
but  it  is  close  enough  so  that  some  improvement  can  be  gained, 

3.5  Application  to  Simulated  Signals 

In  the  case  of  the  generated  received  signal  used  for  the 
example  depicted  in  figures  3-1  through  3-3 i  the  following  parsuneters 
were  used: 

DT  =  5         (microseconds) 

Tq  =  llo5 

T^  =  Tq  +  39.5  ;    nj^  =  8,    A^  =  -0.1 

B^  =  0.8 


When  the  technique  developed  herein  was  applied  to  this  signal,  the 
resulting  estimates  were: 


1 


Si 
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11. '+8785 
8 

-0.096033 
Tq  +  39.5198 
(0.803^17,  -0.000881) 
0.0006864  radians 


=  0.803^17  -  jO. 000881 


The  phase  error  in  this  case,  ^a  ~  i^n  "^  ^0'  ^^   0.039  in 
degrees.  At  a  carrier  frequency  of  100  kHz,  this  corresponds  to  a 
time  error  of  about  1.1  nanoseconds.  The  procedure  was  repeated 
for  a  variety  of  cases  and  the  results  are  presented  in  Table  3-1 • 


^0 

^0 

^1 

^1 

^0 

^1 

^1 

?o  (^^^) 

^„  error 
nanoseconds 

0 

10.0 

40.0 

.8 

10.0004 

39.9996 

(.800009, 
.000004) 

-.000043 

-.068 

0 

10.0 

39.5 

.8 

9.9404 

39.4^04 

(.810026, 
-.00015) 

.000077 

.044 

0 

10.0 

38.0 

.8 

9.8O9I 

37.9863 

(.811575, 
-.00036) 

.000251 

.400 

0 

11.5 

40.0 

.8 

11.5326 

39.9815 

(.7913^, 
-.00017) 

.000125 

.195 

0 

11.5 

39.5 

.8 

11.4878 

39.5198 

(.803417, 
-.00088) 

.000686 

1.095 

0 

11.5 

38.0 

.8 

11.3366 

37.9698 

(.803411, 
,00019) 

-.000227 

-.361 

0 

8.5 

40.0 

.8 

8.5685 

39.9902 

(.795314, 
.00024) 

-.000239 

-.380 

0 

8.5 

39.5 

.8 

8. 5189 

39.5029 

(.799695, 
.00054) 

-.000464 

-.740 

0 

8.5 

38.0 

.8 

8.3650 

37.9386 

(.808735, 
.00013) 

-.001570 

-2.50 

Table  3-1  Results  of  Simulation:  No  Distortion,  One  Echo,  Zero  Phase 

(Note:  For  convenience,  the  T^  column  in  the  tables  actually  represents 
skywave  aurrival  time  relative  to  the  groundwave;  i.e.,  (Table  T. )  =  T^-T„) 


The  results  shown  In  the  table  reflect  the  performance  of  the 
estimators  on  signals  which  were  purely  real  to  begin  with.  The 
phase  errors  would  therefore  define  the  precision  available  in  the 
loop  shown  in  figure  2-3  in  which  the  homomorphic  filter-estimator 
serves  as  the  phase  sensitive  detector. 

The  results  of  the  simulation  for  non-zero  phase  cases  are 
shown  in  Table  3-2.  What  can  be  seen  in  these  cases  is  that  the 
Bj 's  are  essentially  -.19723  radians  out  of  phase  with  the  true  valueSo 
This  can  be  explained  by  examining  the  Fourier  Transform  of  p(n): 

P(eJ")  =  B^  +  B  e-J^l  =  B^  (1  +  ^  e-J^"^l) 


0 


0 


In  other  words,  all  magnitudes  and  phases  in  the  complex  cepstriim 
are  measLired  relative  to  the  groundwave  magnitude  and  phase. 


^0 

Tq 

^1 

^1 

^0 

?o  (^^'i) 

^_  error 
nanosec 

.19723 

8.5 

^0.0 

.8 

8.5687 

39.9900 

(.779802, 
-.155877) 

. 196285 

-1.50 

.19723 

8.5 

39.5 

,8 

8.5169 

39.5030 

(.783872, 
-.157803) 

.197302 

.112 

.19723 

8.5 

38.0 

,8 

8.3671 

37.9392 

(.792381, 
-.160622) 

.202399 

8.08 

Table  3-2  Results  of  Simulation j  No  Distortion,  One  Echo,  Non-zero  Phase 


The  cases  in  which  the  skywaves  are  larger  than  the  groundwave 
are  considered  next.  As  discussed  in  Appendix  B,  exponential  weighting 
must  be  used  in  these  cases.  The  only  thing  that  changes  in  the 

AJ 

estimation  procedure  is  that  ^(n),  the  a  priori  estimate,  must  also 
show  the  effects  of  the  exponential  weighting.  This  is  easily 
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accomplished  as  shown  in  Appendix  A,  The  results  are  presented  in 


Table  3-3.  The  extra  column,  B| ,  gives  the  effective  skywave  coefficient 
relative  to  the  groundwave  after  exponential  weighting. 


^0 

^0 

^1 

^1 

Bl 

^0 

^1 

k 

^0  (^^^) 

^^err 

0 

10.0 

40.0 

(2.0, 
2.0) 

(.51332, 
."51332) 

10.0001 

39.9999 

(.513279, 
.513268) 

.000024 

.039 

0 

10.0 

39.5 

(2.0, 
2.0) 

(.52213. 
.52213) 

9.9419 

39.^99 

(.523293, 
.533561) 

-.00581 

-9.30 

0 

10.0 

38.0 

(2.0, 
2.0) 

(.54945, 
.54945) 

9.8493 

38,0031 

(.547588, 
.560679) 

-.00357 

-5.7 

.78540 

10.0 

40.0 

(2.0, 
2.0) 

.725 

10.0000 

39.9841 

(.719902, 
.001483) 

.784204 

-1.90 

.78540 

10.0 

39.5 

(2.0, 
2.0) 

.74 

9.9450 

39.4468 

(.748461, 
.001661) 

.784109 

-2.05 

.78540 

10.0 

38.0 

(2.0, 
2.0) 

.775 

9.8156 

37.9970 

(.785317, 
.000722) 

.784751 

-1.02 

Table  3-3  Results  of  Simulation:  No  Distortion,  One  Echo,  Skywave 
Larger  than  Groundvfave 


In  the  cases  presented  thus  far,  either  the  groundwave  and  the 
skywave  were  in  phase,  or  one  of  the  two  was  at  zero  phase.  For 
completeness,  Table  3-4  presents  the  results  for  different,  non-zero 
phases o 


^0 

T 
^0 

^1 

^1 

^1 

^0 

A/ 
^1 

?0  (rad)^error 

.19723 

10.0 

40.0 

(i.-i) 

(.39233, 

-.58828) 

10.0002 

39.9993 

i .392375, 
-.58821  ) 

.197271  .064 

.19723 

10.0 

39.5 

(i,-i) 

(.39233, 
-.58828) 

9.9412 

39.4472 

(.388902, 
-.598^+2) 

.201864  7.36 

.19723 

10.0 

38.0 

(i,-!-) 

(.39233 
-.58828) 

9.8I34 

37.9994 

(.38774, 
-.59804) 

.199813 

4.10 

Table  3-4  Results  of  Simulation:  No  Distortion,  One  Echo,  Skyrmve 
and  Groundwave  at  Arbitrary  Phases. 
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In  all  of  the  above  cases,  the  final  estimate  of  the  groundwave 
cepstrum  sequence,  s(n),  was  short  pass  filtered  after  I50  micro- 
seconds before  taking  the  inverse  DFT  to  measure  J^q.  As  discussed 
in  Appendix  A,  this  is  done  to  facilitate  filtering  in  the  multiple 
skywave  case.   Table  3-5  shows  the  results  of  the  estimation  scheme 
when  applied  to  a  signal  with  two  skywaves.  The  results  are  not 
quite  as  good  as  in  the  single  skywave  ca'se.  Nevertheless,  it  is 
felt  that,  if  vjarranted,  more  elegant  estimation  algorithms  can  be 
developed.  The  point  of  the  presentation  thus  fau:  is  to  demonstrate 
that,  in  the  distortion  free,  noiseless  case,  measurement  precision 
on  the  order  of  a  few  nanoseconds  is  possible.  The  effects  of 
signal  distortion  can  now  be  examined 0 
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Chapter  k 
DISTORTED  SIGNAI5 

The  results  of  the  preceeding  chapter  depend  heavily  on  the 
assumption  that  the  received  signal  is  identical  in  form  to  the 
sequence  f(n)  defined  in  equation  (3.1)! 


f(n)  =  's(n)  +  ^(n)  +  ^(n) 


In  practice,  it  can  be  expected  that  g  (n)  =  g(n)  since  this  component 

3, 

is  merely  due  to  the  small  portion  of  the  linear  phase  not  removed 
before  computing  the  cepstrum.  There  remain  however,  two  soiorces 
of  error  since  s  (n)  and  p  (n)  may  not  be  of  the  expected  form.  The 
sensitivity  of  the  estimators  to  these  errors  v:ill  be  examined  in 
this  chapter, 

4,1  Distorted  Basic  Waveform 

In  the  estimation  scheme  of  Chapter  3  it  was  assumed  that  a  good 
a  priori  estimate  of  s  (n)  vias  available.  What  is  meant  by  good  might 
be  that  some  criterion  such  as  equation  (4.1)  is  met. 


(^.1) 


ejn) 

a^  ^ 


s^(n)  -  s(n) 


<  £  for  all  n  such 

that  s  (n)  7^  0 

3> 


where  £  is  some  appropriate  bound.  Upon  close  examination  however,  it 

can  be  seen  that  what  is  of  importance  is  the  effect  of  the  residual 

se'quence,  t"  (n) ,  on  the  sequence  x(n) , 
s 
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x(n)   =  r(n)  -  s(n)   =  s^(n)  +  g^(n)  +  Pg^(n)  -  s(n) 


=  gjn)  +^Jn)  +e^(n) 


Consequently,  what  is  of  concern  is  the  ratio, 


(^,2)         ■/■  7  \ — ;  A'-/  ^   =   /:-  /  \  ,  A — 7~K        5      for  all  n  of 
^   '  gg^(n)  +  p^(n)      g^(n)  +  p^(n) 

Interest 


As  an  example  of  the  distinction,  the  ratio  of  equation  (^.1) 
might  have  a  value  of  about  0.01  for  some  value  of  n.  As  this  stands 
by  itself,  ^(n)  might  therefore  be  considered  a  reasonably  good 
estimate.  However,  if  s  (n)  is  100  times  greater  than  p  (n)  +  g  (n) 
at  this  value  of  n,  the  ratio  of  equation  (^.2)  would  equal  1.  The 
estimate  might  then  not  be  considered  so  good  -  particularly  if  the 
value  of  n  in  question  happened  to  be  near  n^ . 

Conversely,  if  the  ratio  of  equation  (^.1)  is  quite  large  for 
some  quefrency  value  not  considered  in  the  estimation  procedure,  or, 
for  some  quefrency  at  which  p  (n)  +  g  (n)  is  much  larger  than  s  (n) , 

3<  3.  3. 

there  may  be  no  problem.  The  criteria  for  judging  the  appropriateness 
of  the  a  priori  estimate  would  therefore  be  as  in  equation  (^.3)« 

(i+,3)  ^s^"^  ^<  ^a^"^    '      n  =  -1,  -2,  -3,  -^ 


6g(n)  «  $^(n)    ;      n  a;  n^ 


Distortion  of  s  (n)  will  be  introduced  by  the  propagation  path 
and  by  the  receiver  processing.   It  may  be  assumed  that  the  propagation 
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distortion  of  the  groundwave  is  negligible.  Not  so  negligible 
however,  would  be  the  effects  of  receiver  filters.   In  addition  to 
performing  the  low  pass  filtering  in  the  demodulation  scheme  shown 
in  figure  2-2,  the  filters  axe  called  upon  to  suppress  atmospheric 
noise  and  interference  from  other  navigation  or  communication  systems 
operating  in  the  low  frequency  portion  of  the  spectrum.  Since  no 
general  linear  time  Invariant  filter  can  be  used  to  model  the  effects 
of  these  filters  under  all  conditions  of  operation,  what  will  be 
presented  will  be  a  method  of  analyzing  the  technique  which  could 
then  be  applied  to  any  specific  implementation. 

Suppose  the  received  signal  is  processed  by  a  single  pole  low 
pass  filter  (similar  arguments  can  be  applied  for  double  or  triple 
pole  filters).  A  model  of  the  basic  waveform,  i.e.,  the  signal 
without  echoes,  to  be  sampled  would  then  be, 


s^(t)  =  s(t)  *  h(t) 


where 


h(t)   =  e"^^  u_^  (t) 


so  that         s  (n)  =  s(n)  +  h(n) 


a 


and  £g(n)  =  s^(n)  -  s(n)  =  s^(n)  -  s(n)  =  h(n) 

In  this  case,  the  cepstrum  of  the  filter  response  is  easy  to 
compute : 


H(eJ")  =    E  e-^"  e'J™ 
n=0 


1  -  e-(^^<^") 


-1 
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log  H(e-^") 


oo    -6n 


-jvm 


n=l 


n 


for  p  >  0 


and 


^  (n) 


h(n)  = 


-Pn 


;  otherwise 


;  n  >  0 


The  fact  that  the  cepstrum  is  zero  for  negative  quefrencies  will 
be  true  for  any  minimum  phase  filter  -  any  linear  phase  term  being 
removed  in  the  phase  unwrapping  procedure.  Since  the  Aq  estimate  is 
made  from  the  negative  quefrency  values  of  the  received  signal,  it 
can  be  seen  that  the  accuracy  of  the  estimated  sequence  g  (r)  will 

Si 

not  be  affected  by  the  presence  of  the  h(n)  term.   The  values  of 

A 

€  (n)  for  quefrencies  near  n.  however,  will  cause  problems.  The 
values  of  £„(n)  and  p(n)  in  this  region  are  given  in  Table  4-1  for 
comparison.  The  parameters  of  the  signals  in  this  case  are: 

p   =  a  =  2/65 

n^   =  8 

^   =  -0.1 

=  0.8 


DT 


n 

f„(n) 

P(n) 

n^-2 

. 13857 

-.0348 

n,-l 

.11518 

.0792 

"1 

.09770 

.7796 

n^+1 

.08ij-23 

-.0791 

n^+2 

.07351 

.0472 

Table  4-1  Comparison  of  A  Priori 

Estimate  Error  and  p(n) , 


Gleaxly, 
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is  too  laxge  in  this  case  to  permit  estimation 


of  A-t  a-rid  B.  .  The  problem  can  be  overcome  if  the  response  of  the 
filter  can  be  presumed  to  be  known  reasonably  well.  The  a  priori 
estimate  of  s  (n)  can  then  be  modified: 


so  that 


iM 

= 

s(n) 

+ 

A^ 

h(n) 

iM 

= 

^,(n) 

- 

?(n) 

= 

s(n) 

+ 

h(n) 

h(n) 

h(n) 

-  s(n)  - 


h(n) 


In  this  case,  an  estimate  of  the  form 

0       ;  otherwise 


A/ 

h(n)  = 


-Yn 


n 


;  n  >  0 


where  y  ~  P  might  be  used.  The  error  therefore  would  be, 


fjn)  = 


_1_ 

n 


e-pn  _  ^-yn 


n  >  0 


For  Y  -     o"P»  with  the  other  parameters  as  in  the  preceeding  example, 
the  error  values  would  then  be  as  shown  in  Table  ^-2« 
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n 

£%(") 

P(n) 

"1-2 

.01222 

-.03^8 

n,-l 

.01176 

.0792 

"1 

.01129 

.7796 

n^+1 

.01089 

-.0791 

n^+2 

.010^-8 

.0^72 

Table  4-2  Comparison  of  Revised 
A  Priori  Estimate 
Error  and  p(n) 

Within  this  range  of  quefrencies,  £  (n)  <  p  (n; ,  Further 

s      a 

A 

notice  should  he  taken  of  the  fact  that  the  £    (n)  terms  are  of 

s^  ^ 

relatively  constant  value  and  of  constant  sign.  As  developed  in 
Chapter  3i  the  individual  estimates  of  A^  and  B.  will  therefore 
have  nearly  complementary  errors.  Consequently,  the  total  estimates 
of  A-]  a-nd  B^  are  relatively  unaffected  by  the  incorrect  a  priori 
estimate  of  P. 

If  the  correct  value  of  Y  is  used,  i.e.,  Y  =  P,  but  an  incorrect 
filter  form  is  assumed,  the  results  axe   quite  different.  Suppose 


h(t) 


t  e 


-Pt 


u_^(t) 


It  can  be  shown  that  the  resulting  cepstruni  is, 


h(n)  = 


0 

-  P 


-3n 


n 


n  <  0 
n  =  0 
n  >  0 
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Therefore, 


■K(n) 


0 

-Sn 

6  ^ 

n 


" 

t 
1 

n  <  0 

f 

n  =  0 

• 
r 

n  >  0 

This  will  produce  the  same  situation  as  that  depicted  in  Table  4-1. 
In  summary,  it  can  be  said  that  the  a  priori  estimate  of  s  (n) 

3. 

must  be  changed  to  reflect  the  effects  of  the  responses  of  any  filters 
used.  The  estimation  procedure  of  Chapter  3  will  be  relatively 
insensitive  to  errors  in  filter  parameters,  e.g.,  pole  locations, 
but  fairly  sensitive  to  errors  in  estimated  filter  form, 

4.2  Distorted  Echoes 

If  the  received  signal  is  processed  by  a  linear  time  invariant 
filter,  the  effects  will  be  as  described  above.  Such  a  filter  will 
not  change  the  echoed  property  of  the  signal,  i.e.,  p  (n)  vfill  not 
be  changed.  Distortion  of  p  (n)  can  only  be  introduced  by  differences 
in  the  propagation  paths „  The  distortion  can  be  modelled  as  follows 
in  a  single  skywave  case. 


r(t)   =  s^(t)   +  B^  \(^-Tjl)  *  h(t) 


Hen      =  S^ie^-Gjen 


1  + 


B^  e-J^W-'^^-HCeJ^^) 
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log  R(eJ")   =  log  S  (eJ")   +  log  G  (e"^") 

^,,   B,^  e-Jwk(n^+Ai)       . 

+    Z      i-Xf"-'     -^ ^ H^(e-^") 

k=l  ^ 


for 


Bi  H(e-3") 


<   1 


r(n)  =  s^(n)  +  g^(n)  +  p^(n) 
r 


where       p  (n)  =  F 

3< 


^^,   B,^  e-J"^("l^^)   ^   . 

k=l  ^ 


The  resulting  p  (n)  can  be  anticipated  by  recognizing  that  its 

components  are  inverse  Fourier  Transforms  of  products.  Hence,  what 

was  the  first  term  in  the  expansion  is  now  convolved  with 

h(n)  =  h(t)       .  The  second  terra  is  convolved  with  h(n) 

t=n'DT  ^ 

twice  and  the  third  term  is  convolved  with  h(n)  three  times,  etc. 

The  results  of  these  convolutions  are  combined  to  form  p  (n) o  In 

3. 

the  cases  considered  thus  far,  h(t)  has  been  an  impulse  function 

so  that  the  terms  appeared  intact.  If  h(t)  is  only  approxi- 

sin  X 
mately  impulsive,  the  first  term  will  be  only  slightly 

distorted.  Higher  order  terms  however,  will  be  convolved  with 

h(n)  a  number  of  times  and  will  tend  to  spread  out  -  only  loosely 

resembling  functions. 

In  the  Loran-G  case,  h(t)  is  nesirly  impulsive.  Additionally, 

since  the  cepstrum  is  short  pass  filtered,  the  effects  of  the  many- 

sin  X 
fold  convolutions  which  distort  the  high  order  terms  will 
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only  slightly  be  felt.  The  selection  of  an  appropriate  model  for 
h(t)  is  confounded  by  the  fact  that  it  is  so  neaxly  impulsive  that 
it  has  not  been  possible  to  measure  its  form.  Perhaps  the  only 
correct  model  would  be  that  the  skywave  and  groundwave  are  identical 
to  within  \%   along  the  leading  edges  of  the  pulses. 

Although  an  appropriate  form  of  h(t)  is  not  available,  it  is 
desired  to  test  the  sensitivity  of  the  estimators  to  a  possible 
1%   distortion.  What  will  be  used  in  this  case  is  a  groundwave  of 
the  form  t  e    and  a  skywave  of  the  form  t  e   .   If  a  =  2/65 
and  3  =  2/70,  and  if  the  signals  are  scaled  to  achieve  a  value  of 
1  at  their  respective  peaks  (  2/a  and  2/^),  the  skywave  will  have  a 
magnitude  of  about  ^6%   of  the  magnitude  of  the  groundwave  at  t  =  65. 
(with  both  signals  beginning  in  coincidence).  While  the  h(t)  that 


would  produce  such  distortion,  i.e.,   L 


(S-Kl)/(S40)) 


3 


,  is  not 


claimed  to  be  a  particularly  realistic  model  of  what  actually  occurs, 
it  does  serve  as  a  clearly  laxge  overbound  on  the  amount  of  distortion 
to  be  expected. 

The  case  described  above  was  simulated  and  the  results  are 
presented  in  Table  ^-3o  ^ftiile  the  errors  axe,   understandably, 
higher  than  in  the  cases  presented  in  Chapter  3»  they  are  the  results 
of  distortion  which  is  much  more  than  can  be  expected  to  be  introduced 
by  the  propagation  path.  As  long  as  the  propagation  path  filter  is 
nearly  impulsive,  the  effect  on  the  cepstrum  estimators  is  minimal. 

The  results  presented  in  this  chapter  will  be  discussed  further 
in  Chapter  6,  It  can  be  said  however,  that,  with  the  assumptions  made 
in  the  presentation,  signal  distortion  only  slightly  affects  the 
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performance  of  the  homomorphic  filter  as  used  in  this  application. 


fo 

^0 

^1 

^1 

H 

^0 

\ 

^l 

^0  ^^^ 

^0 

error 

.19723 

10.0 

^0.0 

(2., 
-2.) 

(.32977, 

10.0016 

39.875 

(.317^3, 
-.^•7539) 

.16990 

43.5 

.19723 

10.0 

39.5 

(2., 

(.33627, 
-.50^21) 

9.9506 

39.^^662 

(.3^000, 
-.52213) 

.18663 

16.9 

.19723 

10.0 

38.0 

(2.. 
-2.) 

(.35652, 
-.53^59) 

9.8317 

38.0789 

(.3580^+, 
-.55027) 

.18418 

20.8 

Table  4-3  Results  of  Simulation:  Distorted  Echo 
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Chapter  5 
NOISY  SIGNALS 

In  this  chapter  the  effects  of  noise  on  the  cepstnun  estimators 
will  be  examined.  Before  quantitatively  presenting  the  effects  of 
the  noise,  an  attempt  should  he  made  to  develop  a  theory  which 
outlines  the  results  to  be  anticipated.  The  development  begins 
with  an  explanation  of  a  matrix  notation  which  seems  best  suited 
to  a  clear  presentation  of  the  elementary  aspects  of  the  problem. 

If  a  continuous-time  function,  x(t) ,  is  sampled  every  DT  time 
units,  the  resulting  samples  have  values  x(n'DT).  These  samples 
produce  a  sequence  of  values  denoted  x(n)  as  used  in  the  preceeding 
chapters.  The  values  can  also  be  denoted  x  and  used  to  define  a 
vector  x; 


X  =  (xq,x^,x^, . . .x^_^)        =  (x(0),x(1),x(2),...x(N-1)) 


The  sequence  x  has  a  Discrete  Fourier  Transform,  denoted  x,  : 

N-1     _  .i2TTkn 

X,   =    E  X  e~   N 
/Ok       „  n 
n=0 

The  values  of  the  DFT  can  also  be  represented  as  defining  a  vector 


X2 


In  matrix  notatation,  the  DFT  relation  can  be  expressed  as: 
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or, 


(5.1) 


^0 


^1 


^N-1 


X  =  W  X 
/■J 


1 
1 
1 


W 


W 


.   .   .   1 


.  .  .  w 


y(N-l)  „2(N-1) 


N-1 


•  •  •  w 


2(N-1) 


,(N-1)' 


• 

'^o 

^1 

• 

^2 

• 

• 

• 

.  \-l 

where  W  =  exp(-  ■'\r")» 

It  is  desired  to  trace  the  effects  of  the  noise  from  the  time 
domain  to  the  quefrency  domain.  The  first  step  therefore  is  to 
examine  the  behavior  of  the  DFT  of  a  noise  sequence.   To  facilitate 
the  derivation,  the  noise  sequence  to  be  encountered  will  be  considered 
to  consist  of  samples  of  a  zero  mean,  complex,  white  gaussian  random 
process.  The  appropriateness  of  this  model  is  discussed  in  Appendix  D. 
To  proceed,  make  the  following  definition, 

X  =  V  ;      V  is  a  zero  mean,  complex,  white 

gaussian  random  vector 

9 
The  sequence  x  then  has  the  following  properties  : 


E 


[x]  =  E  [Re(x)]   =  E  [lm(2c)]  =  0 


E     Re(x^)   Re(x) 


£ 


[imCx"^ 


)    Im(x) 


^,-  I  . 


^i 
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I  =  N  X  N  identity  matrix 


E 


Re(x  )  Im(x) 


=  0 


E 


X  xr*   =    O      I 


where  *  denotes  complex  conjugate 


In  short,  x  is  a  vector  of  N  complex  random  variables ,  The  vector 
contains  a  total  of  2N  uncorrelated  zero  mean  gaussian  random 
variables  (GRV's),  The  uncorrelated  property,  for  zero  mean  GRV's 
implies  independence.  The  complex  variables,  x.  ,  can  also  be  represented 
as. 


X.   =   X.  . exp 

1     111   ^ 


ARG(x.) 


where  jx. j  is  a  Rayleigh  distributed  random  variable  and  ARG(x.)  is 
a  random  variable  uniformly  distributed  on  (-ttiTt)  o 

One  way  to  arrive  at  the  statistics  of  x  is  to  argue  along 
the  following  lines.  From  equation  (5»l)» 


x,  =  x„  +  vr  x^  +  VJ   x„ 
'-k      0       1       2 


.....W^(^-^)x,_, 


(5.2) 


=   x-  4-  X-  -i-  X-  H-  ...  +X^_^ 


where  x,'  =  W   x, .  The  purpose  of  defining  this  new  sequence  x'  is 
as  follows: 
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•  I  U  T-^  -^ 

xM   =  |W   X  I 


^^' 


Ix'lis  therefore  Rayleigh  distributed  with  statistics  identical  to 
those  of  Ix.  I  .  Furthermore,  ARG(x.')  =  ARG('ri  "  )  +ARG(x.).  Since 
the  phases  are  computed  as  principal  values,  ARG(x.')  is  also  uni- 
formly distributed  on  (-Tr,Tr).  Consequently,  x.'  has  Rayleigh  magnitude 
and  uniform  phase  values  which  are  statistically  identical  to  those 
of  X.  .  It  follows  that  the  real  and  imaginary  parts  of  x.'  must  also 
be  statistically  identical  to  those  of  x. ,  i.e.,  independent  zero 
mean,  identically  distributed  GRV's,  It  is  also  true  that  the 
elements  of  the  x'  vector  are  independent: 


E 


x'  x'-'<- 
m  n 


E 


,,mk    ,T-nk  „ 
m      n 


=  ^^(ra-n) 


E 


x   X* 

.  m  n 


=  a^a 


mn 


By  the  zero  mean  GRV  property,  this  implies  independence.  The  result 
of  this  intermediate  step  can  be  seen  by  re-writing  equation  (5»2) 
as  below. 


(5.2) 


/5k  ~  x^  +  x]^  +  x^  +  „..  +  x'^_^ 


(5.3) 


Re(xj^)  =  Re(x^)  +  Re(xp  +  ...  +  Re(x^_^) 


Ini(xj^)  =  Im(x^)  +  Im(x')  +  ...  +  Im(xj^_^) 
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The  right  sides  of  equation  (5 .3)  contain  a  total  of  2N  independent 
zero  mean  GRV's,  The  two  sums,  i.e.,  R(x,  )  and  Im(;g,  )  ,  aire  therefore 
independent.  Additionally,  they  each  are  sums  of  N  independent  zero 
mean  GRV's  with  variances  ~p~»  It  follows  then  that  Re(x,  )  and 
Im(x,  )  are  independent  zero  mean,  identically  distributed,  GRV's  vfith 
variances  N  — r—  so  that  Ix,  I  is  Rayleigh  and  ARG(;gj^)  is  uniform. 

It  is  also  possible  to  show  that  the  individual  elements  of 
the  x-vector  are  independent.   This  follows  from  the  orthogonality 
of  the  rows  of  the  W-matrix ,  i.e.,  ^,    and  x  become  orthogonal 
random  variables.  Since  they  are  also  zero  mean  and  gaussian,  they 
are  independent.  To  see  this  more  explicitly,  consider: 


E 


.\^* 


m 


Since     E 


=     E 


-0-^^'^ 


X.    +  W 


2k 


,^,...,v/^(^^-i) 


^-1 


X*  +  W-'"x*  +  ...   +  W-'^^N-l)  ^^^^ 


X.    x*l     =     o^  6.  .,   it  follows   that 


E 


~k 


--m 


E 


Ix^l^  +  W^^--)    |x^|2  +   ...+w(^-^)(^-)|x^_j2 


o       N-1  /,  V 


n=0 

For  k  =  m,  the  value  of  the  above  sum  is  clearly  No  .  For  k  /  m,  the 
sum  consists  of  N  terms,  all  with  equal  magnitudes  and  with  phases 
equally  spaced  from  -tt  to  tt.  The  resultant  sum  is  therefore  zero  so 
that, 


E 


^2  \^   yn(k-m) 
n=0 
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=  Na^  6 


km 


Once  again,  together  with  the  zero  mean  gausslan  property,  this  implies 
independence. 

In  summarj'',  a  zero  mean,  complex,  white  gaussian  noise  sequence 
has  a  DFT  which  is  also  a  zero  mean,  complex,  white  gaussian  noise 
sequence.  All  that  changes,  statistically,  is  that  the  variances 
are  changed  from  u  to  Ng  , 

In  computing  a  DFT,  it  is  sometimes  desired  to  add  zeroes  to  the 
end  of  the  sequence.  Since  the  same  DT  is  used,  this  does  nothing  to 
change  the  bandwidth  of  the  DFT,  However,  since  there  axe   more  samples 
of  the  time  sequence,  the  DFT  will  have  more  samples  in  it.  This  has 
the  effect  of  increasing  the  resolution  of  the  DFT,  The  effect  of  this 
operation  on  the  statistics  of  the  noise  DFT  should  be  examined. 

Suppose  that  M  samples  of  the  noise  sequence,  x  ,  ajre  available 
and  that  an  N  point  DFT  is  computed: 
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The  x-vector  will,  in  general,  have  N  non-zero  elements  which, 

as  above,  cam  be  shown  to  have  real  and  imaginary  parts  which  are 

2 
independent,  zero  mean  GRV's  with  variances  M  —^  ,     Hence  the 

magnitudes  of  each  point  are  Rayleigh  distributed  and  the  phases 

are  uniformly  distributed.  The  sequence,  however,  is  no  longer 

white.   In  fact,  the  x-vector  will  span  the  same  M  dimensional 

complex  (or  2M  dimensional  real)  space  as  that  generated  by  the 

original  sequence  x«  The  autocorrelation  function  of  the  DFT  can 

be  computed  in  a  straightforviard  manner.  From  this  it  can  be 

shown  that,  for  example,  if  n/m  =  ^,  points  in  the  DFT  4  samples 

apart  are  Independent,  If  N/M  is  not  an  integer,  all  the  points 

are  correlated,  although  only  loosely  for  seperations  of  more  than 

about  2N/M. 

It  X  -  s_  +  v_,   where  £  is  a  deterministic  vector  and  v  is  the 

noise  vector  as  previously  described,  it  follows  from  the  linearity 

of  the  DFT  operation  that 


and 
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(s-  +  V,  )  (s*  +  V*) 
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s,  s*  +  Na  6, 
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With  this  background,  consider  the  DFT  of  the  sequence 
X  =  s  +  V  which  may  appear  as  in  figvire  5-l» 
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Flgure  5-1  Frequency  Plot  of  Signal  and  Noise 


The  solid  line  in  the  figure  represents  the  magnitude  of  the 
DFT  of  the  deterministic  sequence.  The  dashed  line  represents  some 
convenient  probabalistic  bound  on  the  magnitude  of  the  DFT  of  the 
noise  as  determined  by  its  Eayleigh  distribution.  What  must  be  con- 
sidered is  the  effect  of  the  noise  on  the  phase  of  the  composite 
signal.  For  values  of  k  such  that  Is,  I  >  jv,  I  ,  the  vector  s\im  may 

IK'     I  JK I 

be  as  shown  in  figure  5-2(a), 
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Figure  5-2  Vector  Diagram  of  Signal  and  Noise 
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In  this  case,  the  composite  phase  will  follow  that  of  s,  reason- 
ably closely.  As  an  example,  if  |s  I  =  2|v,  I,  the  phase  error,  the 
difference  between  arg(s,  )  and  arg(s  +  v  )  is  at  most  about  31  •   The 
phase  curve  is  therefore  noisy;  but  the  noise  is  small  -  less  than 
Tr/2  in  magnitude. 

For  values  of  k  such  that  Is,  I  <  jv  L  the  vector  sum  may  be  as 
shown  in  figure  5-2(b).   The  phase  error  in  this  case  can  be  considered 
to  produce  a  more  noisy  curve  -  errors  as  great  as  tt  in  magnitude. 
However,  the  implications  of  this  situation  are  more  significant  in 
view  of  the  fact  that  the  phase  curve  will  now  have  the  structure  of 
the  noise  phase,  i.e.  white.   This  causes  problems  because  it  is 
presumed  in  the  phase  unwrapping  algorithm  that  the  phase  has  a  smooth, 
or  continuous,  structure.   If  ig,  1  >  |v  I  for  all  k,  every  point  in  the 

TT 

DFT  phase  will  be  within  +  -r  of  the  correct  value  so  that  the  measured 
number  of  branchings  of  the  arctangent  function  defining  the  composite 
phase  will  be,  at  most,  one  different  from  the  correct  value.   If 
Is,  j  <  iy,  L  however,  there  can  be  any  number  of  branchings  of  the  phase 
so  that  a  huge  error  can  be  expected  in  the  corrected  phase  curve. 
This  can  be  summarized  by  considering  the  form  of  the  complex 
log  of  the  DFT. 


log 


s,  +  V, 
--k   '-■k 


=  log  fw  +  log 


^k 


-k 
1  +   ^ 


^  . 


-k     kkl 


log  [s^   +  vj   =  log  v^     +     log 
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log  Yv     + 


rk|     l^kl 


If  the  phase  of  v,  is  uniform  and  white,  the  phase  of 


^k 


will  also 


be  uniform  and  white,  when  measured  via  its  principal  value.  Therefore, 

V, 


'k  . 


is  statistically  equivalent  to 


Re(Y^)  +  j  ImCYj^) 


Consequently, 


the  effect  of  the  log  DFT  operation  for  |s,  l»  |v  i  is  to  filter  the 
noise  as  in  figure  5-3(a),   The  Ijg.j  used  is  that  of  the  Loran-G  signal, 


Figure  5-3  LOG  DFT  Filtering  Effect 
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For  1^,  j  <  ly^  if  "the  composite  log  DFT  is  essentially  that  of  y,  , 
Consequently,  the  resultant  log  DFT  can  be  approximated  as: 


log 


s,  +  V, 


^°S5k 


"l  -^ 


log  V,    +  jARG(y,  ) 


'^k 


•»2 


where  H^  and  H^  are  as  shown  in  figure  5-3  (b)  and  (c)  and  i v,  I  and 
ARG(y,  )  are  the  Rayleigh  magnitude  and  uniform  phase  as  previously 
derived.  The  equivalent  operation  is  depicted  in  figure  5-^» 

The  implications  of  this  problem  can  now  be  seen.   The  output  of 
the  phase  unwrapping  block  of  figure  5-^  will  be  grossly  in  error  if 
the  k_  of  figure  5-3  corresponds  to  a  frequency  less  than  the  sampling 
frequency,  1/2DT,   For  a  given  signal  to  noise  level  therefore,  sampling 
faster  than  a  prescribed  ra.te  has  dire  consequences.   It  is,  however, 
desired  to  sample  the  input  signal  as  fast  as  possible,  obtaining  as 
many  samples  of  the  groundwave  as  possible  before  the  first  skywave 

arrives.  The  number  of  samples  so  obtained  will  define  n^  in  the 

sin  X 
complex  cepstriun,  the  quefrency  seperation  between  the  terms 

of  p  (n) ,  If  this  number  is  reasonable,  e.g.,  k-   or  more,  the  estimators 

Si 

can  be  expected  to  perform  well.  If  it  is  less  than  4,  there  will 
be  too  much  overlap  betvieen  the  terms  of  p  (n)  to  obtain  the  desired 

Si 

precision-.  If  a  greater  bandwidth  is  necessary,  some  sort  of 
smoothing  would  have  to  be  applied  to  the  log  magnitude  and  phase 
curves  to  reduce  the  errors  introduced  by  the  false  branchings  due 
to  noise.  In  such  an  operation,  however,  precision  is  lost  in  the 
cepstruffli 


; 

C 

<« 

r-\ 

':3:| 

> 

^ 

^^ 

• 

\^ 

^y 

^•-^ 

4J>   T-) 

^^^^^ 

^ 

i< 

rH 

"■v,,,,^^^ 

>T 

3    !>. 

^■\^^^ 

V—' 

S     rCl 

^\,^^ 

o 

« 

^^--^^           1 

<; 

^..^.^'''''^ 

rf''"*^ 

>^ 

V^^X 

+ 

§ 

rf<— %, 

0 

o 

a 

P4  W 

W? 

rf   nj 

+ 

^£ 

'S' 

^"^ 

g 

£ 

'■<-< 

1 

u 

K 

< 

T-l 

.r 

X 

CM 

y 

> 

. 

I 

r  ^ 

1 

1-H 

-    J 

< 

g 

/ 

<.  1       ■ 

e 

x: 

e 

-p 

X! 

«    -H 

-p 

<D    M 

, 

X   -iH 

iH     rd 

^ 

0)    M 

pi.   tD 

rH     Cd 

M 

13 

« 

^ 

ft   tD 

W) 

m? 

>^ 

§3 

^ 

bO 

U 

o 

o 

rH 

:5| 

' 

E 

> 

+ 

c 

«= 

c 

>o 

-89- 


-90- 

Gonsequently,  a  search  for  a  compromise  along  the  following  lines 
must  be  made, 

1.  For  high  signal  to  noise  ratio,  determine  the  ranges 
of  values  of  sampling  rate  and  arrival  times  in  vfhich 
acceptable  performance  of  the  cepstrum  estimators  is 
achieved, 

2,  Corresponding  to  values  within  these  ranges  of  sampling 
rate  and  arrival  times,  determine  the  lowest  signal 

to  noise  ratio  for  which  the  estimators  will  function. 
The  minimum  SNR  will  determine  how  many  pulses  will  have  to  be 
sampled  to  obtain  the  desired  performance.  Through  considerations 
such  as  signal  stability,  a  means  of  evaluating  the  technique  is 
available. 

An  example  should  now  be  considered  to  see  how  well  the  theory 
developed  thus  far  fits  the  results.  Signal  power  is  typically 
defined  in  terms  of  the  power  at  the  25  microsecond  sampling  point. 
As  it  turns  out,  this  would  be  consistent  with  most  other  definitions 
of  average  power a  Consider  sampling  every  10  microseconds  for  a 
500  microsecond  duration  and  then  taking  a  512  point  DFT  (n/M  x   lO), 

Here,  the  normalized  signal  power  is,  s    «  0.26,  For  each 

2 

component  of  0  db  noise,  the  variance  is  therefore,  a     =0,13.  For 

+40  db  SNR,  o^   =  .000013,  In  the  DFT,  a^   =  50  a^  =  O.OOO65  so  that 

O  =■  . 02^+5. 

From  the  Raylelgh  distribution,  it  can  be  determined  that 

Is,  I  >  iv,  I  with  confidence  levels  ,99  a-nd  ,86  for  log  s,   «  -3,0 
|/vK|   |A-k|  ■"■^  I    "^    '^'kj 

and  -2.6  respectively  ;d.th  v,  as  above.  From  examination  of  a  noise 
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free  DFT  it  can  be  verified  that  these  values  are  achieved  at 
frequencies  of  about  2?  kHz  and  31  kHz,   The  case  was  simulated  and 
the  resulting  magnitude  and  phase  curves  are  presented  in  figure  5-5* 
From  the  phase  curve  it  can  be  seen  that  the  noise  induced  branching 
begins  at  about  36  kHz, 

For  a  +30  db  SNR,  the  same  confidence  levels  correspond  to 
frequencies  of  about  1?  kHz  and  19.5  kHz,  The  results  of  a  simu- 
lation are  shown  in  figure  5-6  where  it  can  be  seen  that  the  first 
branching  occurs  at  about  23  kHz.  For  +20  db,  the  predicted  values 
are  11  kHz  and  13  kHz  while,  as  shovm  in  figure  5-7 i  the  branching 
occurs  at  about  18  kHz,   If  the  magnitudes  of  these  plots  are 
considered,  it  can  be  seen  that  the  frequencies  at  which  branchings 
could  have  occurred  are  nearly  exactly  as  predicted. 

The  cepstrum  for  the  +^0  db  case  was  computed  with  the  following 
signal  parameters:  0^   =  0,  DT  =  20  |isec,  T„  =  20  |j.sec,  T. -T„  =  80  jisec, 
^•i  ~  izii)*     The  resulting  cepstrum  is  shown  in  figure  5-8»  As  can  be 
seen,  the  error  is  small  so  that  the  estimators  can  work  properly.  As 
an  indication  of  this,  the  estimate  of  the  groundwave  phase  was 
-0,00822  radians,  about  a  13  nanosecond  error. 

In  the  +30  db  case,  the  phase  unwrapping  presents  problems  since 
the  25  kHz  point,  i.e.,  the  k  =  N/4  point,  is  above  the  frequency 
at  which  I y,  1  >  Is,  I,  The  cepstrum  is  shown  in  figure  5-9. 

These  examples  are  felt  to  verify  the  theoretical  limitations 
of  the  technique  as  presented  earlier  in  this  chapter,  A  discussion 
of  the  implications  is  included  in  the  following  chapter. 
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Chapter  6 

CONCLUSIONS 

The  results  of  the  preceeding  chapters  will  now  be  summarized. 

The  cepstrum  of  an  undistorted  received  Loran-C  signal  consists 
of  components  vfhich  are  knovm  to  within  a  few  parameters.   If  the 
noise  level  is  such  that  a  sufficient  sampling  rate  can  be  employed, 
the  following  condition  ensues.  While  there  is  overlap  and  interference 
betvfeen  the  components  in  the  cepstrum,  there  are  ranges  of  quefrencies, 
albeit  small,  in  which  particular  components  dominate  the  sum.   In  these 
regions,  parameters  associated  with  the  dominant  component  can  be 
estimated. 

The  performance  of  the  estimators  under  ideal  noise  conditions 
was  demonstrated  in  Chapter  3«  The  technique  estimates  the  groundwave 
phase  as  well  as  could  be  desired  and  has  the  added  advantage  of 
producing  estimates  of  the  various  other  parameters  of  the  signal 
which  could  be  of  interest. 

The  results  of  Chapter  4  show  that  the  signal  distortion  intro- 
duced by  the  propagation  paths,  of  the  amount  to  be  anticipated  in 
the  Loran-C  case,  has  negligible  effect  on  the  estimators.   It  was  also 
shown  that  the  distortion  produced  by  the  receiver  filters  could  be 
accounted  for  so  that  performance  is  maintained  -  if  the  filter 
responses  are  reasonably  well  known. 

The  results  of  Chapter  5»  although  not  conclusive,  give  an 
indication  of  the  eunount  of  processing  time  required  by  the  scheme. 
In  the  +40  db  SNR  case,  a  bandwidth  of  about  36  kHz  was  usable.  A 
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corresponding  sajnpling  rate  of  about  1^  microseconds  would  produce 
values  of  about  3  and  -,5  for  n^  and  A^  respectively  in  the  case  of 
a  35  microsecond  skywave  arrival  time.  Such  a  value  of  n.  is  a 
limit  for  the  type  of  estimators  used  here.  This  is  due  to  the  fact 
that  the  A^  and  B>  estimators  can  not  make  measurements  on  y(n.-2) 
if  n^  =2  since  the  zero  quefrency  point  is  where  all  of  the 
ground wave  magnitude  and  phase  information  is  mapped  to. 

It  has  been  the  purpose  of  this  thesis  to  consider  only  the 
computational  aspects  of  the  application.  Consequently,  a  conclusive 
statement  on  the  feasibility  of  implementing  the  technique  can  not 
be  made.  What  can  be  done  however,  is  to  cite  some  practical 
limitations  of  the  Loran-G  system  and  to  review  the  results  in 
light  of  them. 

The  major  conclusion  of  the  results  is  that  a  signal  to  noise 
ratio  of  about  +  40  db  is  required.   If  each  individual  pulse  was 
at  a  zero  db  power  level,  10,000  pulses  would  have  to  be  averaged  to 
achieve  the  desired  performance  level.  If  all  eight  pulses  of  a 
group  with  a  100,000  microsecond  repitition  interval  were  used, 
this  would  require  an  averaging  time  of  just  over  2  minutes.  Such 
long  times  require  great  signal  stability  and,  of  course,  preclude 
use  in  a  navigational  application,  i.e.,  on  an  aircraft  or  ship. 
The  stability  could  be  achieved  in  a  stationary  application,  as  in 
a  calibration  or  monitoring  receiver,  but  other  effects  must  be 
considered. 

The  eight  transmitted  pulses  are  not  actually  identical. 
There  is  some  degree  of  droop  in  signal  power  since  it  is  difficult 
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for  signal  transmitters  to  recover  from  the  effects  of  a  large  power 
pulse  generated  only  1000  microseconds  earlier.  The  effects  of 
such  droop  can  be  taken  into  account  somewhat  in  the  sampling- 
detection  schemes  presently  employed.  For  a  homomorphic  filter 
receiver, however,  it  may  be  required  that  the  averaging  be  accomplished 
on  the  first  pulse  of  each  group,  thereby  maJcing  the  averaging  time 
16  minutes. 

In  addition,  there  are  the  effects  of  non-linear  filters  used 
to  suppress  the  non-gaussian  aspects  of  the  atmospheric  noise.  These 
and  other  non-linearities  involved  in  actual  processing  would  have  to 
be  taken  into  account,  especially  with  regard  to  whether  they  change 
the  echoed  nature  of  the  signal  drastically. 

Speculation  on  the  feasibility  of  the  technique  therefore  leads 
to  the  following  conclusions.  The  technique  shows  good  results 
when  applied  to  the  most  general  signals  that  can  be  simulated.   In 
view  of  the  long  averaging  times  and  the  more  complex  nature  of  an 
actual  implementation  however,  the  improvement  in  performance  does 
not  seem  marked  enough  to  conclude  that  it  could  be  effectively 
implemented , 
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Appendix  A 
COMPLEX  GEPSTRUM  DERIVATIONS 

This  Appendix  presents  derivations  of  the  complex  cepstra  of  the 
various  signals  considered  in  this  thesis.  The  complex  cepstrum  of 
the  Loran-C  signal,  in  its  most  general  form,  is  somewhat  complicated. 
It  can,  however,  be  arrived  at  in  the  following  systematic  fashion, 

A-1  Basic  V/aveform,  Arrival  Time  Coincident  with  Sampling  Time 

Consider  first  the  basic  waveform,   s(t)  =  t  e    u  .(t),  sampled 
so  that:  s(n)  =  n  e    u  ^ (n) ,  The  first  step  is  to  compute  the 
Fourier  Transform  of  this  sequence.   It  can  be  most  easily  calculated 
by  utilization  of  the  following  property. 

^^'       f (n) F(eJ")  =    E  f (n)  e'^™ 

n=-oo 

then.  1^  F(e^'')     =    S   -jnf (n)  e'^^ 

and,  ^  F(eJ")  =    E   -n^f  (n)  e'J™ 

dw  n=-°o 


The  term  on  the  right  side  of  the  above  equation  is,  by  definition, 
the  Fourier  Transform  of  -n  f (n) ,  It  follows  then  that, 

n2f(n)  ^  -  ^  F(eJ") 

dw 

To  utilize  this  property,  let  f(n)  =  e  "^  u  ^(n).  Then, 


F(eJ")  =    E  e"^  e'^™ 
n=0 


1  _  e"^°^'^J"^ 


-1 
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and, 


dw 


dw' 


1  -  e 


-(a+jw) 


-1 


-(a+jw)  [^_^  g-(a+jw)] 
"  1_  ^-(an-jw)  J  3 


which  is,  therefore,  the  Foiirier  Transform  of  s(n). 

As  described  in  Appendix  B,  the  computation  method  used  in  the 

—  iw 
implementation  will  remove  the  linear  phase  term,  e   ,  thereby  yielding 

a  computed  transform, 


sien 


3^\      _ 


-g  l  +  e-^^-^J^^) 


The  next  step  is  to  compute  the  complex  natural  logarithm  of 


S(e^"). 


log  sCe*^")   =  -a  +  log 

-(a+jw) 


1+  e 


-(a+jw) 


+  3I0S 


Since  a  >0,  and  therefore, 
expansions  are  appropriate: 


1  -  e 
<  1,  the  following  logarithm 


(a+jw) 


log(l+x)  =   Z  (-1)"-^^-^    , 


n=l 

00    X 


n 


n=l 


H<i 


ix  <  1 


Therefore, 


(A.l) 


log  S(e^")  =  -a  +  E 

n=l 

00 
+  3  S 
n=lL 


,1    -an 
(_l)n+l  _e ^-jTm 


-On 
_e g-jwn 


n 


The  final  step  is  to  evaluate  the  inverse  Fourier  Transform  of 
log  S(e  )  and  to  denote  it  as  's(n),  the  complex  cepstnim  of  s(n).  In 
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this  case,  the  form  of  the  expression  in  equation  (A.l)  allows  the 
inverse  transform  to  be  obtained  by  inspection  since  ^(n)  is  the 
coefficient  of  e     in  its  transform. 


(A. 2) 


^(n) 


0 

t 

n  <  0 

-  a 

> 

n  =  0 

3  4-  (-1)"-*-! 

-on 
e 

t 

n  >  0 

Two  generalizations  can  be  made  on  the  results  obtained  thus  far. 
For  any  constant,  B, 

B  s(n) B  S{e^^) >-  log  B  +  log  S(e^^) (log  B)6(n)  +  s(n) 

so  that  any  coefficient  change,  real  or  complex,  in  the  original  signal 
effects  only  the  n  =  0  point  in  the  cepstrum  -  in  an  additive  fashion. 
Additionally,  consider  the  effects  of  sampling  the  continuous- time 
function,   s(t)  =  t  e    u  ^(t),  at  a  rate  (1/dt).  This  yields, 

s(n)  =  (n.DT)^  g-a(n'DT)  ^__^(n.i)T) 
/^^n2  2  -(a-DT)n    /  X 

Thus,  there  is  a  coefficient  change  and  a  virtual  change  in  a  . 
Therefore,  to  compute  the  cepstrum  of  s(n)  for  various  values  of  DT, 
simply  modify  the  expression  in  equation  (Ae2)  by  changing  the  exponent 
from  a  to  (a-DT)  and  adding  the  term  (log  DT)  6(n), 

A-2  Basic  Waveform  Arrival  Time  Between  Sampling  Times 

The  notation  in  this  case  becomes  a  bit  awkward  and  requires  explana- 
tion. The  expression  5(n-x),  for  non-integer  x  will  be  used.  What  is 
implied  is  that  s(n)  ®  6(n-x)  is  a  sampled  version  of  s(t)  *  6(t-x), 


-10^- 


-0.5  <  A<  0.5 


with  this  notation,  consider 

s'(n)  =  s(n)  @  g(n) 

where  g(n)  =  6(n-A)  • 

S'(eJ")  =S(eJ")  G(eJ") 
S(e^   )  will  be  as  derived  in  the  proceeding  section  and 

G(eJ")  =  e-J"^ 
log  Gie^")   =  -jwA 
The  contribution  to  the  cepstrum  of  s'(n)  due  to  g(n)  is, 

TT 


t( 


log     GCe^'^) 


2lT 


S       -jwA  e"^       dw 


-rr 


2Trn 
A 
TTn 

A. 

IT 


[eJ"^"  (jnTT-  1)   +  e-J'""   (jnir  +  l)  ' 
-nrr + 


2j 


sin  nTT 


Ttn  cos  niT 


n 


For  n  =  0, 


g(n)  = 


■n 


0  - 


TT 

n 


cos  nrr 


(.l)n+l       _A_ 
^     ■^  n 


For  n  =  0,  consider  n  to  be  a  continuous  variable  ajid  consider  the 
limit, 


lim 
n—  0 


sin  nrr       -  rm  cos  nrr 


n 


=       lim 
n--0    L 


TT   COS    nTT 

-    TT   COS    nTT 

2 

+  TT  n  sin  nTT 

2n 

=     lira 
n— OL 


tt2 


sin  nTT 


=  0 
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Conseq^uently, 

Since  s' (n)  =  s(n) 
that, 

a'(n)  = 


,  n  =  0 


(_l)n-^l  ^  .  n  ^  0 
g(n),  so  thafs'Cn)  =  s(n)  +'g(n),  it  follows 


(-1) 


n+1 


n 


a 


i=ll 


n-fl^ 


n 


+  ?  +  (-1) e"^" 

n 


,  n  <  0 
,  n  =  0 

,  n  >  0 


A-3  Multipath  Impulse  Train 


The  impulse  train,  denoted  p(n) ,  consists  of 
(a. 3)     p(n)  =  6(n)  +  B^  6(n-n^)  +  ...  +  B^^  ^(^-\) 

To  begin  the  derivation,  consider  the  case  of  m  =  1,  B^j  <  1, 
p(n)  =  6(n)  +  B.  6(n-n.) 

P(e'^")  =    E   p(n)  e-^""    =  1  +  B^  e'^^^ 
n=-«> 

log  P(eJ")  =  log[l  +  B^  e-J^l] 

(B  e"J""l)^ 


The  cepstrura  is  then, 
p(n)   =     F 


CO 

z 

k=l 
-1 


(-1)^ 


since      IB. |  <  1 


log  P(e^") 


1 

2Tr 


"~ 

TT 

s 

-TT 

oo 
I 

k=l 

(-1) 


k+1   ^1 


^-jwkn^ 


,jwn 


dw 


^(n) 


(A.^)  ^(n) 


CO 

k=l 

CO 

k=l 


CO 

E 
k=l 


(-1)     -^      S      ^ dw 

-Tf 


-lOd- 


k    j2Tr(n-kn.) 

k+1  ^1  ^^^   iT(n-knj^) 
k     'rT(n-kn^) 


TT 


-TT 


(-1) 


If  n^  is  an  integer,  equation  (a.^)  becomes 

(-1)^''^  -^     6(n-knp 


kn)     =       E 
k=l 


so  that  the  echo  contributes  an  impulse  train  to  the  cepstrum  (for 

sin  X 


integer  n.).  In  general,  of  course,  the  more  involved 
of  equation  (a.^)  is  encountered. 


behavior 


The  more  general  form  of  equation  (A.3)  calls  for  more  than  one 
echOo  In  the  case  of  m=2, 

p(n)  =  5(n)  +  B^   6(n-n^)  +  B^  eCn-n^). 
In  this  case, 

P(eJ^)  =  1  +  B^e-J""l  +  B^e-J^^^a 


If 


h  r 


<   1,  then 


log  Tie"^'')   =     log 


•jvm. 


1+  B.e"-^""!  +  B.e'^'^Z 


oo 

=   E 
k=l 


(-1) 


1^    ^   .  -2 


-0™1    a.    T5  <=-J™0^k 
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Hils  is  the  Inverse  Fourier  Transform  of  the  cepstruro, 

2  3 

F[$(n)]  =    B^e-J^"'!    -    \   e-J"^"^   ^   !l_  ^-j^Jn^       . 


By  a  straightforward  extension  of  the  previous  results,  it  can  be  seen 

sin  X 
that  the  cepstrum  vd.ll  consist  of  samples  of  functions 

centered  at  quefrencies  n^,  2n.  ,•>•  kn.  and  n  ,  Zn^,    ...  kn^. 

sin  X 
Additionally,  there  are  samples  of  functions  centered  at  all 

comhinatlons  of  multiples  of  n>  and  n„. 

Clearly,  as  the  number  of  echoes  increases,  the  cepstrum  of  the 

.impulse  train  quickly  becomes  complicated  and  difficult  to  filter.  To 

overcome  some  of  the  complications  however,  use  can  be  made  of  the 

form  of  the  signal  information  transfer  between  the  time  and 

quefrency  domains.  As  Schafer  has  shown,  if  the  cepstrum  of  a 

signal  is  short  pass  filtered  with  a  cutoff  quefrency  of,  for  example, 

ra,.,  the  recovered  waveform  is  Identical  to  the  original  waveform  for 

time  values  up  to  m^-..  Additionally,  for  reasonably  large  values  of  m^. 

the  recovered  viaveform  follows  the  original  signal  very  closely  for  some 

time  after  m^^.  Due  to  the  pulse  nature  of  the  Loran  signal,  this  fact 

makes  a  combination  filter  effective. 
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The  cepstrum  can  be  short  pass  filtered  with  a  cutoff  quefrency 
of  about  150  p,sec.  This  vjill  retain  almost  all  of  the  groundwave 
information  but  remove  much  of  the  complicated  impulse  train  contri- 
bution, A  simple  example  can  illustrate  the  effectiveness  of  this 

approach.  Suppose  n.  «  40  (isec,  n_  «  80  |isec,  and  n^~  120  [isec. 

sin  Y 
If  a  cutoff  quefrency  of  I50  p.sec  is  used,  sampled  ^  functions 

centered  at  the  following  quefrencies  are  retained: 

n^.  2n2 

n^^+n^,  2n^+n2,  'i^n^n^ 
Some  of  the  above  quefrencies  aire  above  the  I50  |isec  cutoff,  but  close 
enough  that  their  effects  may  be  felt  in  the  pass  band.  Other  sampled 
— - —  functions  are  centered  at  quefrencies  far  enough  removed  from 
the  cutoff  that  they  may  be  considered  removed  by  the  short  pass 
filter.  Consequently  only  the  quefrencies  listed  above  have  to  be 
considered  in  the  more  elaborate  filter-estimators  of  Chapter  3t 

A-4  Cepstrum  for  Exponentially  Weighted  Signals 

The  multipath  impulse  train  thus  far  has  been  considered  to 
be 

p(n)  =  6(n-nQ)  +  B^  6(n-nj^)  +  B^  6  (n-n^)  +  ..,  +  Bj^  6(n-nj^) 

vd.th  the  following  condition  met: 


(A.5) 


k 

E 
m=l 

B 
m 

<  1 
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This  ensured  that  the  grounduave  was  the  dominant  component  of  the 
composite  signal  and  allowed  the  expansion  of  log  p(e  )  in  the 
manner  presented.  As  will  be  shown  in  Appendix  B,  this  is  part  of  a  set 
of  conditions  sufficient  to  ensure  a  complex  cepstrum  of  the  form 
desired  for  the  methods  of  this  thesis. 

In  the  interesting  Loran-C  cases  this  condition  is  typically  not 
met.  The  problems  can  be  circumvented  however,  by  exponentially 
weighting  the  signal  in  the  manner  to  be  described.  Suppose 

x(t)  =  s(t)  +  B^  s(t-T^) 

u_^(t)  +  B^  (t-T^f   e-^^^-^l^  u_^(t-Ti) 


.2  -at 
=  t  e 


with 


^1 


>  1 


-3 1 
If  the  signal  is  multiplied  by  e   ,  what  results  is 


'e-^^x(t)  =  tV(°^^^)^  u_^(t)  +  B-  (t.T^)2  e-(^-^)(^-Tl)  u_^(t-Tp 


where  B^     =  B^  e"^'^l 


If  3  is  such  that  Ib?  I  <  1,  the  condition  of  equation  (A. 5)  is 
satisfied.  The  resulting  cepstrum  will  be  as  previously  derived  with 
a  replaced  by  (a+3)  and  B^  replaced  by  B> .  The  cepstrum  can  then  be 
computed  and  filtered  in  the  manner  developed  in  this  thesis.  If  it 

is  desired  to  recover  the  original  signal  after  filtering,  the  filtered 

+6t 
signal  can  be  multiplied  by  e     as  a  final  step. 
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The  results  of  the  above  argviment  can  be  generalized  to  the 

case  of  multiple  echoes  to  obtain  the  requirement  that 

If 

E  B  e  ^  m    <   B„  e^^O 
1  I  ml  0 

m=l  '   ' 
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Appendix  B 

PHASE  UNWRAPPING  CONSIDERATIONS 

As  mentioned  in  Chapter  1 ,  the  computation  of  the  complex 
logarithm  of  the  DFT  of  the  signal  presents  problems,  A  computation  on 
a  point  by  point  basis  via  an  arctangent  power  series  expansion  will  only 
provide  the  principle  value  of  the  desired  phase.  What  would  result 
in  the  case  of  a  typical  Loran-C  signal  might  produce  a  phase  plot  as 
in  figure  B-1  (the  dotted  lines). 


-TT 

-2Tr  -1 
-3tt 

-4tt 


Figure  B-1  Uncorrected  Phase  Curve 
It  is  desired  to  remove  the  discontinuities  induced  by  the  branch- 
ing of  the  arctangent  function  for  the  following  reason.  If  x(n)  = 
s(n)  @p(n),  then  it  is  desired  that 


log  Xie^""^   =  log  sCeJ"")  +  log  P(eJ") 
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In  particular,  the  above  equality  must  hold  for  both  the  real  and 
imaginary  parts  of  the  expression.  Thus, 

log  X(eJ")   =  log  S(eJ")   +  log  P(e-^") 
(B.1)       and, 

atg  xCe-^")  =  arg  SCe"^")  +  arg  PCe"^") 

If,  for  example,  there  is  some  value  of  w  for  which  ajrg  S(e'^  ) 
has  a  value  3n/^,   and  arg  T(e"^   )   has  a  value  n/Z,    then  the  right  side 
of  equation  (B,1)  equals  5tt/^.  However,  arg  X(e"^  ),  as  evaluated  by 
the  arctangent  power  series  will  have  a  value  of  -3Tr/^  so  that  the 
equality  is  violated, 

Schafer  has  described  a  method  for  removing  the  discontinuities. 
After  the  complex  log  of  the  DFT  of  the  signal  is  computed,  a  sequence 
with  imaginary  paxt  appearing  as  the  dotted  lines  in  figure  (B-1) 
results.  The  imaginary  part  of  the  sequence  is  then  scanned  to  detect 
jumps  of  almost  Ztt  in  consecutive  values.  For  appropriate  values  of 
k  at  which  these  jumps  occur,  a  correction  sequence  is  generated. 
This  is  represented  by  the  heavy  lines  in  the  figure.  The  coirrection 
sequence  is  added  to  the  originally  computed  phase  sequence  as  shown 
in  figure  B-2. 

As  can  be  seen  from  the  corrected  phase,  there  is  typically  a 
laJTge  linear  phase  term  involved.  Since  this  term  can  easily  dominate 
the  phase  plot,  and  the  resulting  complex  cepstrura,  it  should  be 
removed  if  it  has  no  effect  on  the  results,  Schafer 's  method  is  to 
examine  the  value  of  the  corrected  phase  curve  for  values  of  k  near  n/2. 
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Figure  B-2  Corrected  Phase  Curve 

The  integer  multiple  of  n   that  is  closest  to  this  value  is  called  the 
linear  phase  term.  From  this  value,  the  straight  line  shovm  in  figure 
B-2  can  be  generated.  If  this  is  subtracted  from  the  corrected  phase 
curve,  a  final  version  of  the  phase,  as  in  figure  B-3  results. 


Figure  B-3  Corrected  Phase  Curve  -  Final  Version 
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The  effect  of  this  last  operation  is  to  multiply  the  DFT  of  the 
signal  by  e  "^  0  where  m„  is  the  integer  linear  phase  term  deduced  in 
the  manner  described  above.  This  is  equivalent  to  convolving  the  signal 
with  6(n-ra„)  -  or  shifting  the  sequence  m„  time  samples.  As  a  final 
step  after  the  filtering  is  accomplished  and  the  inverse  operation 
is  performed,  the  recovered  signal  can  be  re-shifted  this  amount  to 
produce  an  output  with  a  proper  time  base. 

The  corrected  complex  logarithm  is  now  continuous  and  periodic, 
and  the  desired  complex  cepstrum  can  be  computed  from  it.  In  effect, 
what  is  presumed  here  is  that  the  branchings  of  the  phase  curve  are 
caused,  in  a  systematic  fashion  throughout  the  complex  log  sequence, 
by  the  same  linear  phase  term.   This  presumption  should  be  examined 
more  closely.  Consider  the  impulse  train  component  of  a  typical 
echoed  signal, 

k 
p(n)  =   Z  B  6(n-n  ) 
1=0        ^ 

The  Fourier  Transform  of  this  sequence  is 


P(e^")  =  E  B.  e"^™i  where  B,   =  Ib. 

1=0  ^  i     1  ^ 


The  following  analysis  will  consider  the  form  of  P(e'^")  for 


various  values  of  the 


B. 
1 


's,  ^.*s,   and  n, 's.  To  begin  the  analysis, 
let  Bq  =  1,  ^„  =  0,  n^  =  0,  and  k  =  1,  Furtliermore ,  let  B^  be  real, 
positive  and  less  than  1  with  n.  =  an  integer.  What  results  in  this 
case  is. 
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0  <  Bj^  <  1 


P(eJ")  =  1  +  B^e-J""l 

n^  =  an  integer 
A  vector  diagram  representation  of  P(e"^  )  is  els  shovm  in  figure  B-4, 


Figure  B-^  Vector  Representation  of  P(e  ) 


As  w  ranges  from  -rr  to  -Ut   the  term  B>  e  "^  1  traces  out  the  circle 
n.  times.  If  B.  is  complex,  i.e.  B.  =  |B.  |  e'^'^1,  0.   f   0,  the  dashed 
line  of  figure  B-4  begins  (for  k  =  O)  at  an  angle  J2f>  with  the  line 
representing  the  vector  1.  Importantly,  as  long  as  B^  <1,  the  branch 
point  at  +rr  of  the  eirctangent  function  defining  the  phase  is  avoided. 
Consequently  the  magnitude  and  phase  of  P(e  )  are  as  plotted  in  figure 
B-5. 

For  B.   <  1,  A,  the  maximum  value  of  the  phase  in  figure  B-5(b), 
is  always  less  than  Tr/2.  If  n.  is  not  an  integer,  the  B^e  '^  1  vector 
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Figure  B-5  Magnitude  and  Phase  of  P(e*^") 

in  figure  B-4  starts  in  line  with  the  1  vector  but  ends  up,  for  w  =^  , 
at  some  arbitrary  angle.  In  terms  of  figure  B-5,  this  means  the  plot 
is  essentially  the  same  -  except  that  the  w  =  +rT  points  are  no  longer 
nulls  or  peaks  of  the  magnitude  curve,  nor  zeroes  of  the  phase  curve. 
For  B^  =  |b^J  e-^^l,  0^  ^  0,  the  entire  plot  of  figure  B-5  is  shifted 
to  the  right  or  left. 

With  these  preliminaries,  more  general  impulse  trains  can  be 
considered.  Suppose 

p(n)  =  6(n-nQ)  +  B^  dCn-n^-n^), 


Consider  first,  the  case  wherein  n-  =  1,  B^  =  a  real,  positive 
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nuraber  less  than  1,  The  magnitude  of  P(e*'  )  will  be  as  in  figure 
B-5(a).  The  phase  however,  will  be  as  in  figure  B-6(a).  When  the 
two  step  correction  procedure  is  employed  as  previously  described,  what 
results  is  presented  in  figure  B-6(c)  and  is  seen  to  be  identical  to  the 
phase  curve  of  figure  B-5(b). 
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Figure  B-6  Frequency  Representation  of  P(e"^*')  With  Linear  Phase 
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A  vector  diagraun  of  P(e  ),  parametric  in  w,  is  shovm  for  this 
case  in  figure  B-7,  Notice  that  the  total  resultant  number  of  crossings 
of  the  branch  line  is  one,  i.e.  n„.  This  is  true  independent  of  the 
orientation  of  the  aixes  in  the  figure,  i,e»  for  arbitrary  ^q. 


/^6 

^Im[p(eJ") 

\e[p(eJ")] 

X 

r^  J 

\                                    y'lncreasing  w 

Figure  B-7  Vector  Representation  of  P(e"^")  with  Linear  Phase 


The  following  summary  can  be  made  of  the  results  thus  fax.  If 


p^(n)  =  6(n)  +  B^  6(n-n^) 


B^  <  1 
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and, 


then, 


PgCn)  =  sCh-Hq)  +  Bj^  6(n-n^-nQ) 


p^Cn)  =  p^(n)  X  6(n-nQ), 


BJ  <1 


Hq  =  an 

integer 


i.e.,  Pp(n)  is  p^ (n)  shifted  n„  samples  to  the  right.  Since  any  linear 
phase  terra  is  removed  in  the  computation,  p^ (n)  and  Pp(n)  become 
identical  signals  beyond  the  complex  log  stage  in  the  processing.  It 
need  merely  be  remembered  to  shift  the  recovered  signal  back  n^ 
samples  as  a  final  step  in  the  second  case. 
Now  consider  the  case  wherein 

p(n)  =  6(n)  +  B^  6(n-n^)         ^1  ""  l^ll  >  •^'' 
The  resulting  phase  cirrve  is  as  shown  in  figure  B-8,  What  can  be 
seen  is  that  the  linear  phase  term  is  due  to  n^ ,  The  ripples  in 
B~8(c)  are  as  in  figure  B-5(b)  but  reversed  in  sense.  This  can  be 
attributed  to  the  fact  that  the  1  vector  in  P(e  )  is  rotating 
counterclockwise  relative  to  the  linear  phase  term  associated  with  n. , 
Further  light  can  be  shown  on  the  effect  by  considering  the  logarithm 
expansion  of  P(e  )  in  this  case. 


log  P(e^")  =  log 


=  log 


1  +  B^  e-J^ 


B^  e-J"^l 


+  log 


1  +  B^-^  e-'J^^l 


-jwn.  +  log  B.   +   E 
■^         ^  k=l 


B. 


-k 


(.l)k+l  _^eJ"^"l 


■•TT  : 
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Figure  B-8     PCe^^   )    Phase,      B.      >  1 
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With  the  linear  phase  term  removed,  the  cepstrum  becomes 


p(n)  =  (log  B.)  6(n)  +    Z 
^  k=l 


-k 


(.l)k-Hl  ^  J  eJw(n4-kn^)  ^^ 

^  '      2nk 

-or 


=  (log  B.)  6(n)  +    Z 

k=lL 


(-1)    _1_  6(n+knp 


for  an  integer  n^.  This  means  that  the  impulse  train  appears  in  the 
negative  quefrency  portion  of  the  cepstrum.  What  is  further  implied 
is  that  the  largest  component  of  the  signal  contributes  to  the  low 
quefrency  portion  of  the  cepstrum  while  the  smaller  components 
produce  the  impulse  trains. 

In  the  Loran-G  case,  it  is  desired  to  recover  the  groundwave 
(or  at  least  to  treat  it  as  the  primary  component) .  Consequently 
exponential  weighting  is  used  to  ensure  that  the  groundwave  is  the 
largest  signal.  This  also  simplifies  the  filtering  process  since,  in 
a  multiple  echo  case,  in  which  a  middle  echo  is  the  largest,  impulses 
would  appear  in  both  the  negative  and  positive  quefrency  regions. 
Additionally,  various  cross  product  terms  in  the  logarithm  expansion 
would  produce  impulses  throughout  the  cepstrum  (e.g.,  at  n^  -  2n. , 
n-3  -  3n.,  etc.) 

From  figure  p,8(b),  it  can  be  seen  how  easily  the  linear  phase 
component  can  obscure  the  important  information  in  the  phase  sequence, 
Clearly,  it  is  iraportajit  that  this  term  be  removed.  However,  in  a 
multiple  echo  case,  the  discontinuities  need  not  be  caused  merely 
by  the  linear  phase  component  of  the  dominant  signal.  In  terms  of 
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the  logarithm  expansion, 

log  P(e^*")   =  log 


E  B.  e-J""l 
i=0  ^ 


the  requirement  can  be  seen  to  be 

k 


(B.2) 


B, 


E  B,  e-J""i 
1=1  ^ 


to  ensure  that  branching  is  due  only  to  the  lineax  phase  term 
associated  with  n„. 

Schafer  states  that,  after  low  pass  filtering  to  avoid  aliasing, 
the  condition, 


(B.3) 


k 

E 
i=i 


B, 


< 


B 


0 


is  sufficient  to  meet  the  requirement  of  equation  (B,2),  This  is 
certainly  consistent  with  everything  presented  thus  far  and  would  be 
.true  if  an  ideal  analog  low  pass  filter  were  employed.   In  such  a 
case,  the  condition  of  equation  (B.3)  could  be  met  via  exponential 
weighting  as  presented  in  Appendix  A,  However,  for  an  analog  filter 
with  a  finite  passband  to  stopband  transition  width,  the  effects  of 
aliasing  should  be  examined  more  closely. 

As  it  develops,  for  transforms  with  ripples  as  in  figure  B-5(a), 
the  effects  of  aliasing  can  cause  large  amounts  of  distortion  in  the 
phase  which  can  not  be  overcome  by  analog  filtering.  To  show  this, 
a  more  in-depth  development  of  the  DFT  must  be  presented.  As  an  example, 
the  DFT  of  a  sequence  is  represented  by  the  dotted  lines  in  figure 
B-9  (the  sequence  is  e   u  ^ (n)  ) ,  The  solid  lines  in  the  figure 
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\^^^  ^^,.^^—  curve  1 

Figure  B-9  Real  and  Imaginary  Part  of  e~  u  > (n)  DFT 


represent  the  continuous  Fourier  Transform  of  the  corresponding 
continuous  time  function  depicted  in  a  manner  appropriate  for  compari- 
son vrith  the  DFT,  i.e.,  with  the  negative  frequencies  shovm  to  the 
left  of  k  =  N. 

It  would,  of  course,  be  desired  that  the  DFT  would  consist  of 
samples  of  the  appropriate  solid  lines,  i.e.,  samples  of  ciurve  1  for 
k  <  N/2  and  samples  of  curve  2  for  k  >  N/2,  As  is  the  case  however, 
the  DFT  is  actually  samples  of  the  sum  of  the  two  curves  for  all 
values  of  k.  Due  to  the  bandwidth  of  the  signal  and  the  sampling 
rate  employed  in  this  case,  it  appears  that  the  aliasing  has  a  negligible 
effect.  The  only  errors  that  appear  remotely  discernible  are  in  the 
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imaginaxy  part  of  the  DFT  for  values  of  k  near  N/2,  and  even  these 
axe  small. 

If  the  DFT  is  shown  via  its  magnitude  and  phase  as  in  figure 
B-10  however,  the  evaluation  changes.  There  is  a  large  error  in 
the  phase  for  a  wide  range  of  values  of  k  in  the  vicinity  of  N/2, 
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Figure  B-lO  Magnitude  and  Phase  of  e   u  .(n)  DFT 


What  can  be  seen  is  that  for  k  <  n/2,  the  composite  phase  is 
in  the  same  quadrant  as  that  of  curve  1.  For  k  >  N/2,  the  composite 
phase  is  in  the  same  quadrant  as  the  phase  of  curve  2,  This  is  due 
to  the  fact  that  the  magnitude  of  the  continuous  Fourier  Transform  of 
e  \x  . (t)  decreases  monotonically  with  frequency.  Consequently  the 
aliased  signal  (curve  i  for  k  >  N/2,  curve  2  for  k  <  n/2)  is  always 
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smaller  than  the  desired  signal  so  that  the  sum  is  dominated  by  the 
correct  signal  and  the  resulting  phase  is  of  the  correct  sign. 

With  this  background,  the  actual  signals  treated  in  this  thesis 
should  be  considered  to  show  what  can  possibly  go  wrong.  Suppose 

x(n)  =  s(n)  (g)  p(n) 
so  that 

X(eJ")  =  S(e^")  P(eJ") 


S(eJ") 


P(e^") 


exp 


arg  SCe-^  )  +  arg  PCe"^  ) 


JW>| 


Figure  B-11  Magnitude  and  Phase  of  Loran-G  Signal  Transform 
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The  magnitude  and  phase  of  the  continuous  Fourier  Transform  of  s(t), 
the  Loran-G  signal,  are  shovm  in  figure  B-11.  What  is  important 
to  notice  is  that,  for  values  of  k  less  than,  but  nearly  equal  to, 
N/2,  curves  1  and  2  are  nearly  180°  out  of  phase  (modulo  2n). 
Now  suppose  that 

p(n)  =  6(n)  +  B^  6(n-n^~ZL,^) 

"^^^^  B^  =  0.8 

n^  =  an  integer 

A^  •=  -.5 
so  that 

P(eJ")  =  1  +  .8  e-J"^"l+^) 
The  vector  plot  of  this  P(e"^")  is  as  shown  in  figure  B-12. 
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Figure  B-12  Vector  Representation  of  P(e"^^) 
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The  rotating  vector  begins,  for  w  =  0,  at  the  position  shovm  as  line  1, 
As  w  increases  to  -Hr,  the  rotating  vector  moves  in  a  clockwise 
direction  to  the  position  shown  as  line  2.  When  w  is  at  a  value  of 
+(lT-6),  the  vector  is  at  the  position  of  line  3  -  where  the  resultant 
magnitude  of  P(e  )  is  0,2e  As  w  decreases  to  a  value  of  -tt,  the 
rotating  vector  moves  in  a  counterclockwise  direction  to  the  position 
of  line  4,  As  w  continues  to  decrease,  the  vector  arrives  at  the 
position  of  line  5  when  w  =  -(17+6)  -  where  the  resultant  magnitude  of 
P(e"^")  is  1.8. 

The  implications  of  all  of  this  is  as  follows.  Let 


Then, 


X  =  X,  +  X_ 
•^     Id. 


where  X  =  DFT 


X.  =  desired 
signal 

X_  =  aliased 
signal 


X(-^_6)  =  X(eJ('"-^))  +  X(eJ(-^-^>) 


S(eJ^^-^b 


PCeJ^"^-^))! 


exp 


"arg  S(eJ^'"-^))  +  arg  He^^'^'^h 


Sie^^-^-^h 


p(3j(-1T-6)) 


exp 


arg  3(6*^^    ^)     +     arg  P(e  ^    ^) 


To  evaluate  this,  the  follovring  values  should  be  used  in  this  case: 


arg 


S(eJ(^-^b   «  -  ^ 


arg 


S(eJ(-^-^^) 


2 
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argP(eJ(^-^)) 


arg  P(eJ(-^-^)) 


He'^'^-^h      -     .2 


It  follovfs  then  that 


X(-|--6)  :,    .2 


S(e 


J(iT-6))  ,-rf 


+  1.8 


22 


S(eJ(-^-^))  eJ  2 


e  *^2^ 


.2  S(eJ(^-'b   -  1.8  S(eJ(-^-'b 


The  desired  value  of  the  phase  is  about  -  -r—  (or  ~p~) «  If  the 

terra  in  the  brackets  is  greater  than  zero,  this  is  what  will  result. 

3Tr 
If  the  term  in  the  brackets  is  negative,  the  phase  will  be  -  -r — h  tt 

= r(or-'^-'rT=  ~     p     ~     ~   ~?)  •  ^^   other  words  a  branching 

will  have  occurred  -  not  caused  by  the  linear  phase  term  of  the 

ground wave • 

For  the  term  in  the  brackets  to  be  positive,  it  is  required  that 


S(eJ('^-^))>  9  S(eJ(-^-^)) 


The  value  of  6  was  such  that  6(n.+A^)  =  'Tt/2  or,  6  =  g/    .  \» 
For  a  Fourier  Transform  of  a  discrete  sequence,  tt  corresponds  to 


N/2  in  the  DFT.  This  can  be  represented  as  a  frequency  w  so  that 
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2(nj^+Ai) 


w 
c 


w 

c 


>  9 


•/  81       s 


If  n^   =  ^0     (in  many  applications  it  can  be  much  larger),   6  sa  -opr. 
What  is  required  therefore  is  thait 

sCeJ^S  "c) 

Since  S(e'^  )|   has  even  symmetry,  what  is  required  is  that  the 
signal  Fourier  Transform  magnitude  decrease  by  a  factor  of  more 
than  9  as  w  goes  from  .9875  w  to  1.0125  w  •  The  Loran-C  signal, 
by  itself,  decreases  by  a  factor  of  about  l.O?  in  this  band.  Consequent- 
ly, an  analog  filter  vfhich  introduces  an  additional  factor  of  about 
9  decrease  is  needed. 

It  can  be  seen  that  if  the  B^  used  in  the  example  were  0,9» 

the  decrease  would  have  had  to  be  by  a  factor  of  z ^r  =19. 

1  -  »y 

If  n^  were  greater  than  40,  the  transition  band  would  be  narrower. 
Additionally,  greater  values  of  6,  i,e.,  6  =  „/   '^.  \  ,  etc., 
will  also  cause  this  problem  so  that  huge  errors  can  accumulate. 

Before  considering  possible  solutions  to  this  problem,  a  few 
examples  of  the  phenomenon  should  be  presented o  A  2048  pt  FFT  was 
computed  for  a  simulated  Loran-G  pulse.  The  complex  logarithm  of 
the  results  are  shown  in  the  following  figures. 

Figure  B-13  shows  aJi  uncorrected  phase  curve  for  0  <  k  <  N/2 
with  n„  =  0,  no  echo.  Figure  B-14  shows  the  same  for  n„  =  4,  Figure 
B-15  shovfs  the  log  magnitude  and  phase  for  n^  =4,  n^  =  20,  B^  =  .8, 
Notice  that  for  values  of  k  near  n/2  (corresponding  to  w  ~  rr)  there  is 
a  small  amount  of  distortion.  This  represents  the  effects  of  aliasing 
as   in  figure  B-10  -  not  enough  to  cause  branchings. 
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Figure  B-16  shows  the  log  magnitude  and  phase  curves  for  n^-.  ■=  4, 
(n.  +A|)  ~   19«5»  Here  the  problem  can  be  seen  (compare  the  last 
portion  of  the  phase  curve  with  figure  B-6.  Figure  B-l?  shows  a  corrected 
phase  curve  for  n„  =  0,  n>  =10,  B.  =0.8,  The  phase  unwrapping  and 
linear  phase  term  removal  have  been  done  properly.  Figure  B-l 8  shows 
the  corrected  phase  for  n^  =  0,  (n.-hA. )  =  9»5»  B.  =  0,8.  The  large 
error  is  apparent, 

A  method  of  overcoming  this  problem  is  discussed  in  Appendix  C, 
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Appendix  C 

TRANSFORM  PHASE  CORRECTION  METHOD 

This  Appendix  develops  the  technique  used  in  this  thesis  to 
avoid  the  phase  correction  problems  presented  in  Appendix  B, 

One  possible  solution  to  the  problems  would  be  to  use  extra 
exponential  weighting.  As  an  example,  if  the  B^  in  figure  B-12  had 
a  value  of  1/3  instead  of  .8,  the  resulting  vector  representation 
would  be  as  shown  in  figure  C-1 . 


Figure  G-1  Vector  Representation  of  P(e  )  With  Extra 
Exponential  V/eighting 


In  this  case,  the  analog  filter  need  only  introduce  a  factor  of  2 
decrease  in  the  transition  band.  However,  in  view  of  the  fact  that 
the  desired  exponential  weighting  to  be  used  would  be  a  complicated 
function  of  signal  parameters,  and  since  the  transition  band  would 
still  be  very  narrow,  an  alternate  solution  might  be  desirable. 
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In  the  cases  considered  in  this  thesis,  the  basic  waveform  is 

2  -Oct 
known  to  be  of  the  form  t  e     which  has  a  phase  value  of  almost 

-  -^  for  large  w.  From  figijre  B-13  it  can  be  verified  that  the 

phase  at  k  «  N/4  «  512  is  essentially  -1.5rr,   In  figure  B-14  it  can 

be  seen  that  the  phase  at  k  «  512  is: 

,arg  X(nA)  «  1.7  -%r  «  -|^ 
«  _|L  _  z,(V2) 


The  scheme  then  becomes:  examine  the  value  of  the  phase  at 
k  «  N/4,  remove  the  -1.5tt  term  which  is  due  to  the  signal  alone,  and, 
determine  the  closest  integer  multiple  of  11/2  in  what  remains.   This 
integer  is  then  n^.  The  phase  can  then  be  corrected  by  bringing 
the  N/4  point  of  the  phase  up  to  zero  by  using  the  appropriate  linear 
phase  correction.  This  is  done  just  as  explained  in  Appendix  B  except 
that  the  straight  line  in  figure  B-2  goes  from  the  k  =  0  point  to  the 
k  =  N/4  point.  The  rest  of  the  DFT  record  is  then  low  pass  filtered 
in  a  particular  fashion  -  it  is  disregarded.  A  new  N/2  point  sequence 
is  generated: 

^new^^^  =  ^old^^^        k  =  0,1,...,(nA)-1 

^new^^)  =  ^old(-f  -^^^   ^=-f (N/2)-l 

The  effect  of  this  operation  is  to  double  the  effective  sampling 

time  (DT) .  A  standard  digital  low  pass  filter  is  not  used  since  zeroes 

jw 
in  the  sequence  log  X(e  )  imply  unity  magnitude  and  zero  phase  in 


the  corresponding  regions  of  the  X(e  )  sequence.  This  would  produce 
an  even  larger  error  in  the  cepstnun. 

The  effectiveness  of  this  scheme  can  be  seen  by  examining 
figure  B-16  and  noticing  that  the  aliasing  does  not  begin  in  the  phase 
until  well  enough  after  k  =  N/^  to  be  safe.  A  comparison  of  the 
resulting  cepstra  for  the  tvio  methods  of  phase  corrections  is  shown 
in  figures  C-2  and  C-3, 
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Figure  C-2  Cepstriim  Computed  With  Conventional  Phase  Unwrapping 
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Cepstrum  for  no  echo: 

DX  «=  5,0     -  This  'becomes 

the  a  priori  estimate 
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Cepstriim  for  one  echo: 
T^  =  10.0,  T.  =  T^  +  ^0.0 
DT  =  5.0,  B^  =  .S'^Bq 
minus  a  priori  estimate 
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Cepstrum  for  one  echo: 
Tq  =  10.0,  T^  =  Tq  +  39.5 

DT  =  5.0,  B^  =  .8  Bq 
minus  a  priori  estimate 
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Cepstrum  for  one  echo; 
Tq  =  11.5,  T^  =  Tq  +  40.0 

DT  =  5.0,  B^  =  .8  Bq 
minus  a  priori  estimate 
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Figure  C-3  Cepstrum  Computed  With  Modified  Phase  Unwrapping 


-14^ 


Appendix  D 
MORE  GENERAL  NOISE  DFT  DERIVATION 

In  this  Appendix,  the  appropriateness  of  the  white  noise  model 
of  Chapter  5  Is  examined  to  see  if  the  results  depend  critically 
upon  it. 

Suppose  the  actual  noise  can  be  modelled  as  being  obtained  in 
the  manner  shown  in  figure  D-1, 


w(t) 

v(t) 

:yC 

v(n) 

h(t) 

DT 

Figure  D-1  Model  of  Non-white  Noise  Sample  Generator 
Here,  w(t)  is  a  zero  mean,  complex,  white  gaussian  noise  process  and 
h(t)  is  some  appropriate  function  which  can  model  the  noise  structure 
and  any  receiver  filters.  The  analysis  can  proceed  as  follows. 


E 


w(t)  w*(t-n) 


O^   6(h) 


E 


w(t) 


=  0 


For  the  purposes  of  this  analysis  it  will  suffice  to  let 


h(t)  =  P  e-^^  u_^(t) 


so  that 


v(t)  =  S     w(y)  h(t-Y)  dY  =  P  /  w(y)  e'^^^'^^  dy 
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It  follows  then  that 


and 


E 


E 


v(t) 


v„  I     =     E 


P  /     E  [  w(y) 
v(n-DT)        =     0 


,-^(t-Y)^Y     =     0 


Additionally, 


E 


v(t)  v*(ji) 


P  t  ^i 

=  3"^  /  /  E 

-oo  -co 


w(X)  w*(y) 


-P(t-X)  -  (m-y)  ,  ^-v 
e  ^   /  g  ''^  '''  dY  dX 


>2  2 


t  pi 


=  ra"/  /  6(X-y) 


^-P(t-m-X-Y) 


dYdX 


For  samples  of  the  function, 


Ry(n,in)  =  E 


V   V* 

n  m 


=  3^  a2  T  T  6(X-Y)  e-^^-^T^-^T-^-Y)  dY  dX 


— oo   — oo 


For  m  >  n  this  becomes 


R^(n,m)   =  ^2  ^  YV(-°T  -^  -°T  -2^)  dX 


-oo 


o2 


-g(n  DT  +  m  DT  -  2X) 
23 


1  n  DT 


_oo 


P  g    -3(ra-n)  DT 


Similarly,  for  m  <  n, 


r>  (        \  3  O         -3(n-m).DT 


so  that 


R  (1)  .  ^     e-ei-M 


where  1  =  m  -  n 
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The  results  of  Chapter  5  depend  heavily  on  an  analysis  of  the 
variance  of  points  in  the  DFT  -  from  this,  the  noise  magnitude  bound 
is  derived.  Therefore  it  is  important  to  consider  R  (k,k)  «=  E  y,  v* 


E 


V   V* 


E 


r   ^  ,r^    ^  w2k    ^     ^  „k(N-l) 

rO      1       V  +  ...  +  W  '^   '   V  ^^ 


vg  +  W~^  V*  +  „..  +  w"^^'^"^^  v*_^ 


=  N  R^(0)  +  R^(l)  • 

+  R  (2)  • 
v^  >" 


w~^  +  /  +  w"^  +  ...  +  W^ 


w"^  +  w^  +  ...  +  w^ 


=  NR^(O)  +  2  I      (N-m).R^(m)  -cos^ 

ni=l 


D-1 


R^(k.k)  =  ^   N  +  2  ''E\N-m).e-P-^T  •">„-,   ^ 
'^  m=l 


cos 


2TTkm 


The  evaluation  of  the  sum  in  this  expression  can  be  simplified  by 
making  the  approximation  (N-m)  rj  N.  To  see  the  validity  of  this 
approximation,  the  following  would  have  to  be  specified. 


P  «  1.5  X  10- 


DT  =  10  [isec 


p-DT  «  1.5 


and   e"^  °'^  «  .223 


The  h(t)  filter  with  such  a  pole  location  would  have  a  half  power  point 


-11+8- 
at  a  frequency  of  about  25  kHz  so  that  3  has  a  reasoriable  value.  The 
sum  in  equation  (D-l)  then  becomes: 

2   (N-ra)  a  cos  — —     ;      a  ~  ,223 
m=l 

A  plot  of  the  three  terms  in  the  sum  is  shown  in  figure  D-2  for 
N  =  256.  Since  the  rapid  decay  of  the  a  term  will  dominate  the  sum, 
it  can  be  seen  that  (N-m)  ~  N  is  a  clearly  justified  simplification. 
The  sum  to  be  considered  is  then 


„   „   m     2'nK:m 
N   S  a  cos  — z7- 

m=l         ^' 

which  can  be  evaluated  as  in  Jolley 

^   m      ^     1  -  a  cos  9  +  a    cos  (N-l)9  -  a  cos 
E  a  cos  m6  =  ^ ^-r 

m=0  1  -  2a  cos  0  +  a 

For  N  =  256  and  a  «  .223>  a  ^  a.         rj  0,  so  that  the  sum  is 


N-1 
^   m 
L     a  cos  m 

m=l 


N-1 
i-     a  cos  m( 
m=0 


-  1 


a  cos  9  -  a _  a  [cos  9  -  aj 

1  -  2a  cos  9  +  a^     |l  -  a  e""^®|^ 


The  result  is  plotted  in  figure  D-3  for  9  =  — 
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/              2'rTk           2s 
(a  cos  -j^-     -  a  ) 
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---^      1   -  a  e          N 

\                                     k 
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--.223                            ^^^----^.--'^ 

Figure  D-3 
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Figure  D-4  Non-V/hite  Noise  DFT  Variance 
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What  can  be  noticed  is  that  the  zero  crossing  is  at 


or  at 


-1 
G  =  cos   a 


N     -1 
k  =  •^  cos   .223 


.214  N 


(since  e  =  ^  ) 


It  follows  from  equation  (D-1)  that 


R^Ck.k) 


^  L^ 


N-1 


m 


1  +  2  E  a  cos  me 
m=l 


for  a  «  ,223 
2lrtc 


9  = 


N 


This  is  plotted  in  figure  D-4. 

What  can  be  seen  is  that  as  P-DT  increases  (a  decreases) j  the 
plot  flattens  out,  approaching  a  value  of  NR  (O)  as  in  the  white 
noise  case.  The  crossover  point,  i.e.,  the  value  of  k  beyond  which 
the  noise  variance  is  less  than  in  the  vfhite  noise  case,  was  derived 
from, 


N      -1 


-P'DT 


As  P-DT--CO,  e"^'^'^ -0,  so  that  k^  approaches  |-  (-y-)  =  nA. 

As  described  in  Appendix  G,  it  is  the  values  of  the  DFT  for  k  <   n/4 
that  are  to  be  retained  to  compute  the  cepstnim.  Additionally,  the 
sampling  rate  for  any  given  SNR  is  to  be  determined  by  the  noise 
statistics  in  the  DFT  for  values  of  k  neaar  N/4.  What  the  analysis 
has  shown  therefore  is  that  for  k«  N/4,  the  variance  of  the  noise  is 
nearly  the  same  in  both  the  white  and  non-white  cases  so  that  the  model 
used  in  Chapter  5  is  appropriate. 
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